Load packages

x <-
  c(
    "metafor", "kableExtra", "Matrix", "dplyr", "wesanderson", "weightr", "ggplot2", "MuMIn", "RColorBrewer", "tidyverse", "patchwork", "devtools", "R.rsp", "orchaRd")

lapply(x, function(y) {
  # check if installed, if not install
  if (!y %in% installed.packages()[, "Package"])
    install.packages(y)
  
  # load package
  try(require(y, character.only = T), silent = T)
})

#remotes::install_github("itchyshin/orchard_plot", subdir = "orchaRd", force = TRUE, build_vignettes = TRUE)

Read in data

dat <- read.csv("../data/MA_dat-21-06-2023.csv", header = T, sep = ",")

# subset data that have a "variation type". We'll parse based on variation type below
#dat <- dat[-which(is.na(dat$Variation.Type) == TRUE),]
dat <- subset(dat, select = -c(Data_source,Comment)) 
colnames(dat)
##  [1] "ID"                        "Study"                    
##  [3] "Experiment"                "Entered.by"               
##  [5] "Host"                      "Host.type"                
##  [7] "Parasite"                  "Parasite.type"            
##  [9] "Study.type"                "X_gradient"               
## [11] "Gradient.category"         "Trait"                    
## [13] "Trait.type"                "Trait.category"           
## [15] "Gradient.as.predicted"     "Statistics"               
## [17] "Level_C"                   "Level_X"                  
## [19] "N_C"                       "N_X"                      
## [21] "Mean_C"                    "Variation_C"              
## [23] "Lwr_C"                     "Upr_C"                    
## [25] "Mean_X"                    "Variation_X"              
## [27] "Lwr_X"                     "Upr_X"                    
## [29] "Success_C"                 "Success_X"                
## [31] "Variation.Type"            "OR"                       
## [33] "OR_SE"                     "OR_N"                     
## [35] "F_X"                       "F_df"                     
## [37] "Chi2"                      "Chi2_df"                  
## [39] "Slope"                     "R"                        
## [41] "R_N"                       "R_df"                     
## [43] "Slope_SE"                  "Z_X"                      
## [45] "Z_N"                       "Z_df"                     
## [47] "Z_SE"                      "Environment"              
## [49] "Environmental_stage"       "Vector_transmitted"       
## [51] "Final_host"                "Require_intermediate_host"
## [53] "DOI"

Data cleaning

Gradient categories

# Toxin is a type of pollution
dat$Gradient.category <- gsub(pattern = "Toxin", replacement = "Pollution", x = dat$Gradient.category)

# exclude ambiguous/equivocal categories
dat <- dat[which(dat$Gradient.category != "Ecological"),]

# Check gradient categories 
levels(as.factor(dat$Gradient.category))
## [1] "Environment" "Pollution"   "Resource"
# Check pollution experiments are actually chemical pollution
levels(as.factor(dat$X_gradient[dat$Gradient.category == "Pollution"]))
##  [1] "17-estradiol (E2)"                           
##  [2] "Alexandrium tamarense toxin"                 
##  [3] "alpha-cypermethrin"                          
##  [4] "Atrazine"                                    
##  [5] "Atrazine Herbicide"                          
##  [6] "Cadmiun"                                     
##  [7] "Carbaryl"                                    
##  [8] "Carbaryl Pesticide "                         
##  [9] "Copper contamination"                        
## [10] "Deltamethrin"                                
## [11] "diatomaceous earth"                          
## [12] "Diatomaceous earth"                          
## [13] "Diatomaceous Earth (DE)"                     
## [14] "Fipronil"                                    
## [15] "Fipronil Insecticide"                        
## [16] "Flupyradifurone"                             
## [17] "Funguscides"                                 
## [18] "Glyphosate"                                  
## [19] "Glyphosate exposure"                         
## [20] "Imidacloprid"                                
## [21] "Insecticide"                                 
## [22] "Lead Exposure"                               
## [23] "Malathion Insecticide Exposure "             
## [24] "Naphthenic acids"                            
## [25] "Noctural light"                              
## [26] "Noise"                                       
## [27] "Pesticide"                                   
## [28] "pesticide lambda (k)-cyhalothrin"            
## [29] "Pesticides (thiacloprid and tau-fluvalinate)"
## [30] "Thiacloprid"                                 
## [31] "Thiacloprid Insecticide"                     
## [32] "Thiamethoxam"                                
## [33] "Wastewater treatment plant effluents"
# Nocturnal light and noise should be "Environment"
dat$Gradient.category[dat$X_gradient == "Noctural light"] <- "Environment"
dat$Gradient.category[dat$X_gradient == "Noise"] <- "Environment"

# Check Environment experiments are not chemical pollution
levels(as.factor(dat$X_gradient[dat$Gradient.category == "Environment"]))
##  [1] "Artificial light at night"  "Diurnal Temperature Range "
##  [3] "Drought"                    "Environmental enrichment"  
##  [5] "Humidity"                   "Humidity and Temperature " 
##  [7] "Hypoxia"                    "Light exposure"            
##  [9] "Noctural light"             "Noise"                     
## [11] "Oxygen"                     "pH"                        
## [13] "Plant species diversity"    "Salinity"                  
## [15] "Seasonality (Temperature) " "Temperature"               
## [17] "tidal stress"               "Water Temperature"
# Check resource limitation experiments are actual resources
levels(as.factor(dat$X_gradient[dat$Gradient.category == "Resource"]))
##  [1] "Food"                         "Food "                       
##  [3] "Food availability"            "Food level"                  
##  [5] "Low food"                     "Nitrogen and phosphorus conc"
##  [7] "Nutrional stress"             "Pollen limitation"           
##  [9] "Protein (yeast) limitation"   "Selenium supplementation"    
## [11] "Starvation"                   "Starvation time"             
## [13] "Sugar limitation"

Transmission mode

Following R1 suggestion to discuss and/or include as moderator.

for (i in 1:nrow(dat)) {
  if (is.na(dat$Environmental_stage[i]) == F &
      is.na(dat$Vector_transmitted[i]) == F  &
      is.na(dat$Require_intermediate_host[i]) == F) {
    if (dat$Environmental_stage[i] == "Present" |
        dat$Vector_transmitted[i] == "Yes" |
        dat$Require_intermediate_host[i] == "Yes") {
      dat$Tansmission[i] <- "Indirect"
    } else {
      dat$Tansmission[i] <- "Direct"
    }
  } else {
    dat$Tansmission[i] <- NA
  }
}

Median, min, max and n to mean and SD

non_par <- dat[which(dat$Variation.Type == "Median_Min_Max"),]

# Estimate mean for control and treatment from median, maximum and minimum
non_par$Estimated_Mean_C <- (non_par$Lwr_C + 2 * non_par$Mean_C + non_par$Upr_C) / 4
non_par$Estimated_Mean_X <- (non_par$Lwr_X + 2 * non_par$Mean_X + non_par$Upr_X) / 4

# Estimate SD from range and n
# First, read in Xi_N table from Wang et al. (2014) BMC Medical Research Methodology
Xi_table <- read.csv("../scripts/Xi_N_Table.csv", header = T, sep = ",") 

# Estimate SD using the approximation for each N under 50
for (i in 1:length(non_par$ID)) {
  non_par$Estimated_SD_C[i] <-
    (non_par$Upr_C[i] - non_par$Lwr_C[i]) / Xi_table$Xi_N[which(Xi_table$N == non_par$N_C[i])]
  non_par$Estimated_SD_X[i] <-
    (non_par$Upr_X[i] - non_par$Lwr_X[i]) / Xi_table$Xi_N[which(Xi_table$N == non_par$N_X[i])]
}

# overwrite dat
dat$Mean_C[which(dat$Variation.Type == "Median_Min_Max")] <- non_par$Estimated_Mean_C 
dat$Mean_X[which(dat$Variation.Type == "Median_Min_Max")] <- non_par$Estimated_Mean_X

dat$Variation_C[which(dat$Variation.Type == "Median_Min_Max")] <- non_par$Estimated_SD_C
dat$Variation_X[which(dat$Variation.Type == "Median_Min_Max")] <- non_par$Estimated_SD_X

dat$Variation.Type[which(dat$Variation.Type == "Median_Min_Max")] <- "SD"

# Sanity check
levels(as.factor(dat$Variation.Type))
##  [1] ""                    "CI"                  "Cloppler-Pearson CI"
##  [4] "GLMM"                "GM_SD"               "IQ"                 
##  [7] "logistic95CI"        "Median_IQ"           "OR"                 
## [10] "SD"                  "SE"                  "t"                  
## [13] "Wald CI"             "Z"                   "Z_reg"

Median and IQ range to mean

non_par_2 <- dat[which(dat$Variation.Type == "Median_IQ"),]

# Estimate mean for control and treatment from median and IQ range
non_par_2$Estimated_Mean_C <- (non_par_2$Lwr_C + non_par_2$Mean_C + non_par_2$Upr_C) / 3
non_par_2$Estimated_Mean_X <- (non_par_2$Lwr_X + non_par_2$Mean_X + non_par_2$Upr_X) / 3

# overwrite dat
dat$Mean_C[which(dat$Variation.Type == "Median_IQ")] <- non_par_2$Estimated_Mean_C 
dat$Mean_X[which(dat$Variation.Type == "Median_IQ")] <- non_par_2$Estimated_Mean_X

dat$Variation.Type[which(dat$Variation.Type == "Median_IQ")] <- "IQ"

# Sanity check
levels(as.factor(dat$Variation.Type))
##  [1] ""                    "CI"                  "Cloppler-Pearson CI"
##  [4] "GLMM"                "GM_SD"               "IQ"                 
##  [7] "logistic95CI"        "OR"                  "SD"                 
## [10] "SE"                  "t"                   "Wald CI"            
## [13] "Z"                   "Z_reg"

Subset effects with workable “variation types”

nrow(dat)
## [1] 1194
accepted_var <- c("SE", "SD", "CI", "Wald CI", "IQ", "OR")
dat <- dat[which(dat$Variation.Type %in% accepted_var),]

nrow(dat)
## [1] 1144
# store sample size in an object
N <- nrow(dat)

Getting effect sizes

Calculate OR and approximate variance from contingency tables

for (i in 1:N) {
  # no correction need if none of the categories are 0
  if (dat$Variation.Type[i] == "OR" &
      dat$Success_C[i] > 0 &
      dat$Success_X[i] > 0 &
      (dat$N_C[i] - dat$Success_C[i]) > 0 &
      (dat$N_X[i] - dat$Success_X[i]) > 0) {
    dat$OR[i] <-
      (dat$Success_X[i] * (dat$N_C[i] - dat$Success_C[i])) / (dat$Success_C[i] * (dat$N_X[i] - dat$Success_X[i]))
    dat$Log.OR[i] <- log(dat$OR[i])
    # approximate variance
    dat$Log.OR.v[i] <-
      1 / dat$Success_X[i] + 1 / (dat$N_X[i] - dat$Success_X[i]) + 1 / dat$Success_C[i] + 1 / (dat$N_C[i] - dat$Success_C[i])
  }
  else {
    # if at least one category is 0, we apply Yate's correction to avoid dividing by 0
    if (dat$Variation.Type[i] == "OR" &
        (dat$Success_C[i] == 0 |
        dat$Success_X[i] == 0 |
        (dat$N_C[i] - dat$Success_C[i]) == 0 |
        (dat$N_X[i] - dat$Success_X[i]) == 0)) {
      dat$OR[i] <-
        ((dat$Success_X[i] + 0.5) * (dat$N_C[i] - dat$Success_C[i] + 0.5)) / ((dat$Success_C[i] + 0.5) * (dat$N_X[i] - dat$Success_X[i] + 0.5))
      dat$Log.OR[i] <- log(dat$OR[i])
      # approximate variance
      dat$Log.OR.v[i] <-
        1 / (dat$Success_X[i] + 0.5)  + 1 / (dat$N_X[i] - dat$Success_X[i] + 0.5) + 1 / (dat$Success_C[i] + 0.5) + 1 / (dat$N_C[i] - dat$Success_C[i] + 0.5)
    }
    else{
      # if there is no OR data, make these columns NA
      dat$OR[i] <- NA
      dat$Log.OR[i] <- NA
      dat$Log.OR.v[i] <- NA
    }
  }
}

Calculate SD for comparisons with continuous normally distributed variables

# create SD variables
dat$SD_C <- vector(length = N)
dat$SD_X <- vector(length = N)

for (i in 1:N) {
  # Calculate SD if SE is reported
  if (dat$Variation.Type[i] == "SE") {
    dat$SD_C[i] <- dat$Variation_C[i] * sqrt(dat$N_C[i])
    dat$SD_X[i] <- dat$Variation_X[i] * sqrt(dat$N_X[i])
  } else {
    # calculate SD if a 95% confidence interval for a normal distribution is reported
    # this also applies to the Wald confidence interval for proportions, as it is a normal approximation to the binomial
    if (dat$Variation.Type[i] == "CI" |
        dat$Variation.Type[i] == "Wald CI") {
      # calculate SE from lower and upper limits
      temp_C1 <- (dat$Upr_C[i] - dat$Mean_C[i]) / 1.96
      temp_C2 <- (dat$Mean_C[i] - dat$Lwr_C[i]) / 1.96
      
      temp_X1 <- (dat$Upr_X[i] - dat$Mean_X[i]) / 1.96
      temp_X2 <- (dat$Mean_X[i] - dat$Lwr_X[i]) / 1.96
      
      # average digitized/recorded values of SE (because we have two) and transform to SD
      dat$SD_C[i] <-
        mean(abs(c(temp_C1, temp_C2))) * sqrt(dat$N_C[i])
      dat$SD_X[i] <-
        mean(abs(c(temp_C1, temp_C2))) * sqrt(dat$N_X[i])
    } else {
      # approximate SD if IQ range is reported
      if (dat$Variation.Type[i] == "IQ") {
        dat$SD_C[i] <- (dat$Upr_C[i] - dat$Lwr_C[i]) / 1.35
        dat$SD_X[i] <- (dat$Upr_X[i] - dat$Lwr_X[i]) / 1.35
      }
      else {
        # if SD is reported leave as such
        if (dat$Variation.Type[i] == "SD") {
          dat$SD_C[i] <- dat$Variation_C[i]
          dat$SD_X[i] <- dat$Variation_X[i]
        } else {
          # if there is no appropriate data, make these columns NA
          dat$SD_C[i] <- NA
          dat$SD_X[i] <- NA
        }
      }
    }
  }
}

Get effect sizes and variances

# within groups standard deviation
dat$S_within <-
  sqrt(((dat$N_C - 1) * dat$SD_C ^ 2 + (dat$N_X - 1) * dat$SD_X ^ 2) / (dat$N_C + dat$N_X - 2))

# standardized effect size
for (i in 1:N) {
  # standardized mean difference
  # we can't include data without variance
  if (is.na(dat$S_within[i]) == FALSE & dat$S_within[i] > 0) {
    dat$d[i] <- (dat$Mean_X[i] - dat$Mean_C[i]) / dat$S_within[i]
    dat$var.d[i] <-
      (dat$N_C[i] + dat$N_X[i]) / (dat$N_C[i] * dat$N_X[i]) + (dat$d[i] ^ 2) /
      (2 * (dat$N_C[i] + dat$N_X[i]))
  }
  else{
    # if only F statistic is reported - no cases of this yet!
    if (dat$Variation.Type[i] == "F_X") {
      dat$d[i] <-
        sqrt(dat$F_X[i] * (dat$N_X[i] + dat$N_C[i]) / (dat$N_C[i] * dat$N_X[i]))
      dat$var.d[i] <-
        (dat$N_C[i] + dat$N_X[i]) / (dat$N_C[i] * dat$N_X[i]) + (dat$d[i] ^ 2) /
        (2 * (dat$N_C[i] + dat$N_X[i]))
      
    } else{
      # if Z and N are reported for a regression analysis
      if (dat$Variation.Type[i] == "Z_reg") {
        r <- dat$Z_X[i] / sqrt(dat$Z_N[i])
        dat$d[i] <- (2 * r) / (sqrt(1 - r ^ 2))
        Vr <- (1 - r ^ 2) ^ 2 / (dat$Z_N[i] - 1)
        dat$var.d[i] <- (4 * Vr) / (1 - r ^ 2) ^ 3
      }
      else{
        # if t and N are reported for a regression analysis
        if (dat$Variation.Type[i] == "t") {
          r <- sqrt((dat$Z_X[i] ^ 2) / (dat$Z_X[i] ^ 2 + dat$Z_df[i]))
          dat$d[i] <- (2 * r) / (sqrt(1 - r ^ 2))
          Vr <- (1 - r ^ 2) ^ 2 / (dat$R_N[i] - 1)
          dat$var.d[i] <- (4 * Vr) / (1 - r ^ 2) ^ 3
        }
        else{
          # If Z is reported from a comparison of two independent groups
          if (dat$Variation.Type[i] == "Z") {
            dat$d[i] <-
              sqrt(abs(dat$Z_X[i]) * sqrt(dat$N_C[i] + dat$N_X[i]) / (1 - sqrt(
                dat$Z_X[i] ^ 2 * (dat$N_C[i] + dat$N_X[i]) ^ -1
              )))
            dat$var.d[i] <-
              (dat$N_C[i] + dat$N_X[i]) / (dat$N_C[i] * dat$N_X[i]) + (dat$d[i] ^ 2) /
              (2 * (dat$N_C[i] + dat$N_X[i]))
          }
          else{
            # converting log OR
            if (is.na(dat$Log.OR.v[i]) == FALSE &
                dat$Log.OR.v[i] > 0) {
              dat$d[i] <- dat$Log.OR[i] * (sqrt(3) / pi)
              dat$var.d[i] <- dat$Log.OR.v[i] * 3 / pi ^ 2
            } else {
              # for now, leave other types of effects as NA
              dat$d[i] <- NA
              dat$var.d[i] <- NA
            }
          }
        }
      }
    }
  }
}

Correct for sample size and get g

# sample size correction factor
for (i in 1:N) {
  # if odds ratio or comparison between means
  if (is.na(dat$N_C[i]) == F) {
    dat$J[i] <- 1 - 3 / (4 * (dat$N_C[i] + dat$N_X[i] - 2) - 1)
  } else {
    # if correlation
    if (is.na(dat$R_N[i]) == F) {
      dat$J[i] <- 1 - 3 / (4 * (dat$R_N[i] - 2) - 1)
    }else {
      dat$J[i] <- NA
    }
  }
}
# corrected effect size
dat$g <- dat$J * dat$d

# and variance
dat$var.g <- (dat$J ^ 2) * dat$var.d

Flip effect sizes for remaining mortality effects

mort <- grep(pattern = "[M,m]ort", x = dat$Trait)

for(i in mort){
  dat$g[i] <- -dat$g[i]
}

Prepare data

Filter experiments with both epidemiological and infected demographic effects

# For now, get rid of NAs
dat <- dat[(is.na(dat$d) == F),]

# How many experiments do we have?
Exps <- unique(dat$Experiment)

# A data frame of experiments that included Infected demographic and Uninfected demographic effects
dat_1 <- data.frame()

# A data frame of experiments that included Infected demographic and Epidemiological effects
dat_2 <- data.frame()

# populate data sets with corresponding experiments
for (i in Exps) {
  tmp <- dat[which(dat$Experiment == i), ]
  if ("Uninfected demographic" %in% tmp$Trait.category &
      "Infected demographic" %in% tmp$Trait.category) {
    dat_1 <- rbind(dat_1, tmp)
  }
  if ("Epidemiological" %in% tmp$Trait.category &
      "Infected demographic" %in% tmp$Trait.category) {
    dat_2 <- rbind(dat_2, tmp)
  }
}

dat_1 <- dat_1[which(dat_1$Trait.category == "Uninfected demographic" | dat_1$Trait.category == "Infected demographic"),]

dat_2 <- dat_2[which(dat_2$Trait.category == "Epidemiological" | dat_2$Trait.category == "Infected demographic"),]

dat_2$Trait.type <- factor(dat_2$Trait.type, levels = c("Prevalence", "Intensity", "Fecundity", "Survivorship"), ordered = T)

Variance-covariance of sampling errors matrix (for multiple comparisons)

This code produces a variance-covariance of sampling errors matrix, n x n, with n = number of effect sizes. It requires the dataset (dat), Hedge’s g (g), and the variance of Hedge’s g (var.g) to be defined above. It uses “Experiment”, “Level_C”, “Trait”, “N_C”, and “N_X” columns.

Identify negative eigenvalues

First, we find out which studies have negative eigenvalues in their variance-covar matrix. This means covariance between treatment effects with the same control is larger than variance and it’s an indicator of suspiciously low variance and potentially an incorrectly reported N or an error while digitizing data.

For data set 1

studies_1 <- unique(dat_1$Study)
studies_to_check <- c()
negative_eigen <- c()

for (k in studies_1) {
  dat_study <- dat_1[which(dat_1$Study == k),]
  varcovmat = matrix(0,
                     nrow = dim(dat_study)[1],
                     ncol = dim(dat_study)[1])
  
  for (i in 1:dim(dat_study)[1]) {
    for (j in 1:dim(dat_study)[1]) {
      if (i == j) {
        varcovmat[i, j] = dat_study$var.g[i]
      } else{
        if (dat_study[i, "Experiment"] == dat_study[j, "Experiment"] &
            dat_study[i, "Level_C"] == dat_study[j, "Level_C"] &
            dat_study[i, "Trait"] == dat_study[j, "Trait"] &
            dat_study[i, "Trait.category"] == dat_study[j, "Trait.category"]) {
          varcovmat[i, j] = 1 / dat_study[i, "N_C"] + dat_study$g[i] * dat_study$g[j] /
            (dat_study[i, "N_C"] + dat_study[i, "N_X"] + dat_study[j, "N_X"])
        }
      }
    }
  }
  
  val <- eigen(varcovmat)
  for (m in 1:length(val$values)) {
    if (val$values[m] < 0) {
    studies_to_check <- append(studies_to_check, k)
    negative_eigen <- append(negative_eigen, val$values[m])
    }
  }
}
    
Neg_eigen_1 <- data.frame(study = studies_to_check, eigen = negative_eigen)

Neg_eigen_1 %>%
  kbl() %>%
  kable_material(c("striped", "hover"), full_width = F)
study eigen
Shostak et al. 2015. Journal of Parasitology -1.0728652
Ashraf et al. 2017 Environ Sci Pollut Res -0.0311074
Ashraf et al. 2017 Environ Sci Pollut Res -0.0459004
Ashraf et al. 2017 Environ Sci Pollut Res -0.0506383
Ashraf et al. 2017 Environ Sci Pollut Res -0.1843072
Ashraf et al. 2017 Environ Sci Pollut Res -0.2369469
Ashraf et al. 2017 Environ Sci Pollut Res -0.3778444
Ashraf et al. 2017 Environ Sci Pollut Res -0.5191928
Ashraf et al. 2017 Environ Sci Pollut Res -0.5801272
Ashraf et al. 2017 Environ Sci Pollut Res -0.6831878
Ashraf et al. 2017 Environ Sci Pollut Res -0.7934325
Ashraf et al. 2017 Environ Sci Pollut Res -0.8845345
Ashraf et al. 2017 Environ Sci Pollut Res -0.9065495
Ashraf et al. 2017 Environ Sci Pollut Res -0.9289586
Ashraf et al. 2017 Environ Sci Pollut Res -0.9316950
Ashraf et al. 2017 Environ Sci Pollut Res -0.9469193
Ashraf et al. 2017 Environ Sci Pollut Res -1.0721254
Ashraf et al. 2017 Environ Sci Pollut Res -1.1539142
Ashraf et al. 2017 Environ Sci Pollut Res -1.2820301
Ashraf et al. 2017 Environ Sci Pollut Res -1.3782808
Ashraf et al. 2017 Environ Sci Pollut Res -1.4589050
Schluter-Vorberg et al. 2017. Aquatic Toxicology -0.0496863
Schluter-Vorberg et al. 2017. Aquatic Toxicology -0.0753315
Waki & Yoshinaga 2015. FishPathology -0.1246226
Waki&Yoshinaga. 2013. Fish Sci -0.0097697
Waki&Yoshinaga. 2013. Fish Sci -0.0487760
Waki&Yoshinaga. 2013. Fish Sci -0.1770843
Pochini and Hoverman. 2017. Environmental Pollution -0.1229350
Pochini and Hoverman. 2017. Environmental Pollution -1.1461176
Civitello et al. 2020. Proc B -0.5611703
Civitello et al. 2020. Proc B -3.8834534

For data set 2

studies_2 <- unique(dat_2$Study)
studies_to_check <- c()
negative_eigen <- c()

for (k in studies_2) {
  dat_study <- dat_2[which(dat_2$Study == k),]
  varcovmat = matrix(0,
                     nrow = dim(dat_study)[1],
                     ncol = dim(dat_study)[1])
  
  for (i in 1:dim(dat_study)[1]) {
    for (j in 1:dim(dat_study)[1]) {
      if (i == j) {
        varcovmat[i, j] = dat_study$var.g[i]
      } else{
        if (dat_study[i, "Experiment"] == dat_study[j, "Experiment"] &
            dat_study[i, "Level_C"] == dat_study[j, "Level_C"] &
            dat_study[i, "Trait"] == dat_study[j, "Trait"] &
            dat_study[i, "Trait.category"] == dat_study[j, "Trait.category"]) {
          varcovmat[i, j] = 1 / dat_study[i, "N_C"] + dat_study$g[i] * dat_study$g[j] /
            (dat_study[i, "N_C"] + dat_study[i, "N_X"] + dat_study[j, "N_X"])
        }
      }
    }
  }
  
  val <- eigen(varcovmat)
  for (m in 1:length(val$values)) {
    if (val$values[m] < 0) {
    studies_to_check <- append(studies_to_check, k)
    negative_eigen <- append(negative_eigen, val$values[m])
    }
  }
}
    
Neg_eigen_2 <- data.frame(study = studies_to_check, eigen = negative_eigen)

Neg_eigen_2 %>%
  kbl() %>%
  kable_material(c("striped", "hover"), full_width = F)
study eigen
Bradley et al. 2019. Plos One -0.0015336
Bradley et al. 2019. Plos One -0.0250523
Bradley et al. 2019. Plos One -0.0379281
Bradley et al. 2019. Plos One -0.1081902
Souto et al. 2015. VetMicrobiol. -0.5585385
Carrington et al. 2013. PLOS NegTropDiseases -0.0050917
Carrington et al. 2013. PLOS NegTropDiseases -0.0603175
Carrington et al. 2013. PLOS NegTropDiseases -0.1012185
Ashraf et al. 2017 Environ Sci Pollut Res -0.0311074
Ashraf et al. 2017 Environ Sci Pollut Res -0.0506383
Ashraf et al. 2017 Environ Sci Pollut Res -0.2369469
Ashraf et al. 2017 Environ Sci Pollut Res -0.2917917
Ashraf et al. 2017 Environ Sci Pollut Res -0.3005591
Ashraf et al. 2017 Environ Sci Pollut Res -0.3576572
Ashraf et al. 2017 Environ Sci Pollut Res -0.3844600
Ashraf et al. 2017 Environ Sci Pollut Res -0.3890634
Ashraf et al. 2017 Environ Sci Pollut Res -0.3998072
Ashraf et al. 2017 Environ Sci Pollut Res -0.5191928
Ashraf et al. 2017 Environ Sci Pollut Res -0.5801272
Ashraf et al. 2017 Environ Sci Pollut Res -0.6662958
Ashraf et al. 2017 Environ Sci Pollut Res -0.6831878
Ashraf et al. 2017 Environ Sci Pollut Res -0.7429232
Ashraf et al. 2017 Environ Sci Pollut Res -0.7522409
Ashraf et al. 2017 Environ Sci Pollut Res -0.7660915
Ashraf et al. 2017 Environ Sci Pollut Res -0.7758519
Ashraf et al. 2017 Environ Sci Pollut Res -0.7934325
Ashraf et al. 2017 Environ Sci Pollut Res -0.8008515
Ashraf et al. 2017 Environ Sci Pollut Res -0.9216851
Ashraf et al. 2017 Environ Sci Pollut Res -0.9233991
Ashraf et al. 2017 Environ Sci Pollut Res -0.9253113
Ashraf et al. 2017 Environ Sci Pollut Res -0.9316950
Ashraf et al. 2017 Environ Sci Pollut Res -0.9398650
Ashraf et al. 2017 Environ Sci Pollut Res -1.0283485
Ashraf et al. 2017 Environ Sci Pollut Res -1.0694983
Ashraf et al. 2017 Environ Sci Pollut Res -1.0721254
Ashraf et al. 2017 Environ Sci Pollut Res -1.1094691
Ashraf et al. 2017 Environ Sci Pollut Res -1.2178041
Ashraf et al. 2017 Environ Sci Pollut Res -1.3472667
Ashraf et al. 2017 Environ Sci Pollut Res -1.3815538
Ashraf et al. 2017 Environ Sci Pollut Res -1.4370017
Ashraf et al. 2017 Environ Sci Pollut Res -1.5357111
Ashraf et al. 2017 Environ Sci Pollut Res -1.5532399
Ashraf et al. 2017 Environ Sci Pollut Res -1.7623978
Ashraf et al. 2017 Environ Sci Pollut Res -1.9073762
Ashraf et al. 2017 Environ Sci Pollut Res -2.3818336
Ashraf et al. 2017 Environ Sci Pollut Res -2.6483026
Ashraf et al. 2017 Environ Sci Pollut Res -6.8202809
Ashraf et al. 2017 Environ Sci Pollut Res -7.0495394
Ashraf et al. 2017 Environ Sci Pollut Res -8.3833406
Schluter-Vorberg et al. 2017. Aquatic Toxicology -0.0496863
Schluter-Vorberg et al. 2017. Aquatic Toxicology -0.0753315
Waki&Yoshinaga. 2013. Fish Sci -0.0077688
Waki&Yoshinaga. 2013. Fish Sci -0.0487760
Waki&Yoshinaga. 2013. Fish Sci -0.1770843
Purcell et al. 2016 J Fish Diseases -0.1184091

Based on these eigenvalues, we looked back at the studies/experiments and determined if they had issues with reporting N or variance. We decided to remove two such studies for data set 1 (Ashraf et al. 2017 and Shostak et al. 2015) and one such study for data set 2 (Ashraf et al. 2017).

As we found no errors or causes of concern for other studies with barely negative eigenvalues, we are forcing the var-covar matrices to the nearest positive definite values, while keeping the diagonals (var) to the estimated values and modifying only the expected covariances.

Also for now we are excluding Experiment 121 from Civitello et al 2020 Proc. B. This study has a relatively low variance for the control of the uninfected demographic fecundity effects, which results in a negative eigenvalue.

# exclude suspicious studies with highly negative eigenvalues
dat_1 <- dat_1[-grep("Ashraf et al. 2017 Environ Sci Pollut Res", dat_1$Study),]
dat_1 <- dat_1[-grep("Shostak et al. 2015. Journal of Parasitology", dat_1$Study),]
dat_1 <- dat_1[-grep(121, dat_1$Experiment),]

varcovmat_1 = matrix(0, nrow = dim(dat_1)[1], ncol = dim(dat_1)[1])

for (i in 1:dim(dat_1)[1]) {
  for (j in 1:dim(dat_1)[1]) {
    if (i == j) {varcovmat_1[i,j] = dat_1$var.g[i]}else{
      if (dat_1[i, "Experiment"] == dat_1[j, "Experiment"] & dat_1[i, "Level_C"] == dat_1[j, "Level_C"] & dat_1[i, "Trait"] == dat_1[j, "Trait"] & dat_1[i, "Trait.category"] == dat_1[j, "Trait.category"]) {
        varcovmat_1[i,j] = 1/dat_1[i,"N_C"] + dat_1$g[i]*dat_1$g[j]/(dat_1[i,"N_C"] + dat_1[i,"N_X"] + dat_1[j,"N_X"]) 
      }
    }
  }
}

# correct negative eigens in the few studies with large covar to var ratios
varcovmat_1_PD <- nearPD(varcovmat_1,  keepDiag = TRUE)


# exclude suspicious studies with highly negative eigenvalues
dat_2 <- dat_2[-grep("Ashraf et al. 2017 Environ Sci Pollut Res", dat_2$Study),]

# Carrington et al. 2013 has extremely large sampling variance for a parasite prevalence effect. DOUBLE CHECK!!!
dat_2 <- dat_2[-grep("Carrington et al. 2013. PLOS NegTropDiseases", dat_2$Study),]

varcovmat_2 = matrix(0, nrow = dim(dat_2)[1], ncol = dim(dat_2)[1])

for (i in 1:dim(dat_2)[1]) {
  for (j in 1:dim(dat_2)[1]) {
    if (i == j) {varcovmat_2[i,j] = dat_2$var.g[i]}else{
      if (dat_2[i, "Experiment"] == dat_2[j, "Experiment"] & dat_2[i, "Level_C"] == dat_2[j, "Level_C"] & dat_2[i, "Trait"] == dat_2[j, "Trait"] & dat_2[i, "Trait.category"] == dat_2[j, "Trait.category"]) {
        varcovmat_2[i,j] = 1/dat_2[i,"N_C"] + dat_2$g[i]*dat_2$g[j]/(dat_2[i,"N_C"] + dat_2[i,"N_X"] + dat_2[j,"N_X"]) 
      }
    }
  }
}  

# correct eigenvalues in the few studies with large covar to var ratios
varcovmat_2_PD <- nearPD(varcovmat_2, keepDiag = TRUE)

Create a new column with a category of vertebrate or invertebrate (keeping this for now)

dat_1 <- mutate(dat_1, Host.type.2 = case_when(
  Host.type == "Fish" ~ "Vertebrate",
  Host.type == "Arthropod" ~ "Invertebrate",
  Host.type == "Amphibian" ~ "Vertebrate",
  Host.type == "Mollusc" ~ "Invertebrate"))

dat_2 <- mutate(dat_2, Host.type.2 = case_when( 
  Host.type == "Fish" ~ "Vertebrate",
  Host.type == "Arthropod" ~ "Invertebrate",
  Host.type == "Amphibian" ~ "Vertebrate",
  Host.type == "Mollusc" ~ "Invertebrate",
  Host.type == "Reptile" ~ "Vertebrate",
  Host.type == "Bird" ~ "Vertebrate",
  Host.type == "Mammal" ~ "Vertebrate"))

Numbers mentioned in the methods

The Methods section describes various features of the data, indicating the number of studies/effects with those features. Numbers are calculated here.

# First, get all the curated data
final_dat <- merge(dat_2, dat_1, all = T )
nrow(final_dat)
## [1] 891
# do we have any observational studies left after filtering?
levels(as.factor(final_dat$Study.type))
## [1] "Experimental"
# Number of studies that report median instead of mean
c(unique(non_par$Study), unique(non_par_2$Study)) %in% final_dat$Study
## [1] TRUE TRUE TRUE TRUE
length(unique(non_par$Study)) + length(unique(non_par_2$Study)) 
## [1] 4
length(non_par$ID) + length(non_par_2$ID) 
## [1] 13
# number of studies with only data range
length(non_par$ID)
## [1] 8
length(unique(non_par$Study))
## [1] 1
# number of studies with only IQ
length(non_par_2$ID)
## [1] 5
length(unique(non_par_2$Study))
## [1] 3
# number of studies with OR
length(unique(final_dat$Study[final_dat$Variation.Type == "OR"]))
## [1] 66
# number or studies w/ multiple treatment levels or measuring time points
k = c()
for(i in 2:nrow(final_dat)){
 if(final_dat$Trait[i] == final_dat$Trait[i-1] &
    final_dat$X_gradient[i] == final_dat$X_gradient[i-1] &
    final_dat$Host[i] == final_dat$Host[i-1] &
    final_dat$Parasite[i] == final_dat$Parasite[i-1]) {
   k <- append(k, final_dat$Experiment[i])
    }
}
length(unique(k))
## [1] 110
# studies with negative eigen values
sum(unique(c(Neg_eigen_1$eigen,Neg_eigen_2$study)) %in% final_dat$Study == T)
## [1] 5
# how many experiments
length(unique(final_dat$Experiment))
## [1] 122
# how many studies include more than one experiment
k = c()
for(i in 2:nrow(final_dat)){
 if(final_dat$Experiment[i] != final_dat$Experiment[i-1] &
    final_dat$Study[i] == final_dat$Study[i-1]) {
   k <- append(k, final_dat$Study[i])
    }
}
length(unique(k))
## [1] 21
# how many host species per taxonomic groups
Host_count <- final_dat %>% group_by(Host.type) %>%
  summarise(N_taxa = length(unique(Host)))

Host_count %>%
  kbl() %>%
  kable_material(c("striped", "hover"), full_width = F)
Host.type N_taxa
Amphibian 21
Arthropod 20
Bird 2
Fish 13
Mammal 1
Mollusc 13
Reptile 1
# how many parasite taxa per taxonomic groups
Par_count <- final_dat %>% group_by(Parasite.type) %>%
  summarise(N_taxa = length(unique(Parasite)))

Par_count %>%
  kbl() %>%
  kable_material(c("striped", "hover"), full_width = F)
Parasite.type N_taxa
Bacteria 14
Fungus 6
Helminth 12
Myxozoa 1
Protozoan 8
Virus 23

Analysis

Question 1 - Is there a differencial response to stressors between infected and uninfected host?

Model selection

# First, we evaluate some code that generates helper functions needed so that metafor and MuMIn could interact as necessary
eval(metafor:::.MuMIn)

# Now, we take the full model and fit the rest of the models and examine those models whose AICc value is no more than 10 units away from that of the best model
full_model.1 <- rma.mv(
    g,
    V = varcovmat_1_PD$mat,
    mods = ~ Trait.category*Trait.type*Gradient.category, 
    random = list( ~ 1 | Experiment,  ~ 1 | ID),
    data = dat_1,
    method = "ML"
  )

model_selection.1 <- dredge(full_model.1, trace = 2) 
## Fixed term is "(Intercept)"
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=                                                                     |   1%
  |                                                                            
  |=                                                                     |   2%
  |                                                                            
  |==                                                                    |   2%
  |                                                                            
  |==                                                                    |   3%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |====                                                                  |   5%
  |                                                                            
  |======                                                                |   9%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |===========                                                           |  16%
  |                                                                            
  |=============                                                         |  18%
  |                                                                            
  |=================                                                     |  24%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |==========================                                            |  37%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |==================================                                    |  49%
  |                                                                            
  |===================================================================== |  99%
subset(model_selection.1, delta <= 10, recalc.weights = FALSE)
## Global model call: rma.mv(yi = g, V = varcovmat_1_PD$mat, mods = ~Trait.category * 
##     Trait.type * Gradient.category, random = list(~1 | Experiment, 
##     ~1 | ID), data = dat_1, method = "ML")
## ---
## Model selection table 
##     (Int) Grd.ctg Trt.ctg Trt.typ Grd.ctg:Trt.ctg Grd.ctg:Trt.typ
## 22      +       +               +                               +
## 24      +       +       +       +                               +
## 56      +       +       +       +                               +
## 1       +                                                        
## 32      +       +       +       +               +               +
## 5       +                       +                                
## 3       +               +                                        
## 7       +               +       +                                
## 64      +       +       +       +               +               +
## 2       +       +                                                
## 6       +       +               +                                
## 4       +       +       +                                        
## 39      +               +       +                                
## 8       +       +       +       +                                
## 128     +       +       +       +               +               +
##     Trt.ctg:Trt.typ Grd.ctg:Trt.ctg:Trt.typ df   logLik   AICc delta weight
## 22                                           8 -555.340 1127.1  0.00  0.405
## 24                                           9 -554.802 1128.1  1.02  0.243
## 56                +                         10 -554.791 1130.2  3.11  0.086
## 1                                            3 -562.492 1131.0  3.98  0.055
## 32                                          11 -554.509 1131.7  4.66  0.039
## 5                                            4 -561.842 1131.8  4.72  0.038
## 3                                            4 -561.889 1131.9  4.82  0.036
## 7                                            5 -561.316 1132.8  5.73  0.023
## 64                +                         12 -554.457 1133.8  6.69  0.014
## 2                                            5 -561.802 1133.8  6.70  0.014
## 6                                            6 -561.026 1134.3  7.21  0.011
## 4                                            6 -561.168 1134.6  7.49  0.010
## 39                +                          6 -561.313 1134.8  7.78  0.008
## 8                                            7 -560.480 1135.3  8.19  0.007
## 128               +                       + 14 -553.535 1136.2  9.14  0.004
## Models ranked by AICc(x)
# relative importance values for the predictors can be obtained with:
sw(model_selection.1)
##                      Trait.type Gradient.category Gradient.category:Trait.type
## Sum of weights:       0.88       0.84              0.79                       
## N containing models:    14         14                 6                       
##                      Trait.category Trait.category:Trait.type
## Sum of weights:       0.48           0.12                    
## N containing models:    14              6                    
##                      Gradient.category:Trait.category
## Sum of weights:       0.06                           
## N containing models:     6                           
##                      Gradient.category:Trait.category:Trait.type
## Sum of weights:      <0.01                                      
## N containing models:     1

Full Model

We fit the full model to show there is no effect of infection status

Q1.f <-
  rma.mv(
    g ~ Trait.category:Trait.type:Gradient.category -1,
    V = varcovmat_1_PD$mat,
    random = list(~ 1 | Experiment,  ~ 1 | ID),
    data = dat_1,
    method = "REML"
  )

pal <- RColorBrewer::brewer.pal(9,"Spectral")

Q1.f_OP <-
  orchard_plot(
    Q1.f,
    xlab = "Standardised mean difference (Hedge's g)",
    transfm = "none",
    angle = 0,
    alpha = 0.7
  ) +
  scale_fill_manual(values = rep(pal, each = 2)) + scale_colour_manual(values = rep(pal, each = 2)) +
  theme(legend.position = c(0.25, 0.01), axis.title = element_text(size = 15), axis.text.y = element_text(size = 15)) +
  scale_y_discrete(
    labels = c(
      "Infected:EE:fecund",
      "Uninfected:EE:fecund",
      "Infected:EE:surv",
      "Uninfected:EE:surv",
      "Infected:CP:fecund",
      "Uninfected:CP:fecund",
      "Infected:CP:surv",
      "Uninfected:CP:surv",
      "Infected:RL:fecund",
      "Uninfected:RL:fecund",
      "Infected:RL:surv",
      "Uninfected:RL:surv"
    )
  )

Q1.f_OP

Best Model: Effects in survival vs fecundity under different types of stress, combining infected and uninfected

Q1 <-
  rma.mv(
    g ~ Trait.type:Gradient.category -1 ,
    V = varcovmat_1_PD$mat,
    random = list( ~ 1 | Experiment, ~ 1 | ID),
    data = dat_1,
    method = "REML"
  )

summary(Q1)
## 
## Multivariate Meta-Analysis Model (k = 384; method: REML)
## 
##    logLik   Deviance        AIC        BIC       AICc  ​ 
## -548.7330  1097.4659  1113.4659  1144.9451  1113.8562   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed      factor 
## sigma^2.1  0.2560  0.5060     78     no  Experiment 
## sigma^2.2  0.3607  0.6005    384     no          ID 
## 
## Test for Residual Heterogeneity:
## QE(df = 378) = 217601.2265, p-val < .0001
## 
## Test of Moderators (coefficients 1:6):
## QM(df = 6) = 46.7815, p-val < .0001
## 
## Model Results:
## 
##                                                      estimate      se     zval​ 
## Trait.typeFecundity:Gradient.categoryEnvironment      -0.3255  0.2429  -1.3404 
## Trait.typeSurvivorship:Gradient.categoryEnvironment   -0.5161  0.1016  -5.0776 
## Trait.typeFecundity:Gradient.categoryPollution        -0.1406  0.3153  -0.4458 
## Trait.typeSurvivorship:Gradient.categoryPollution     -0.3510  0.1477  -2.3767 
## Trait.typeFecundity:Gradient.categoryResource         -0.8637  0.2161  -3.9974 
## Trait.typeSurvivorship:Gradient.categoryResource       0.0025  0.2050   0.0121 
##                                                        pval    ci.lb    ci.ub 
## Trait.typeFecundity:Gradient.categoryEnvironment     0.1801  -0.8015   0.1505 
## Trait.typeSurvivorship:Gradient.categoryEnvironment  <.0001  -0.7153  -0.3169 
## Trait.typeFecundity:Gradient.categoryPollution       0.6557  -0.7586   0.4775 
## Trait.typeSurvivorship:Gradient.categoryPollution    0.0175  -0.6404  -0.0615 
## Trait.typeFecundity:Gradient.categoryResource        <.0001  -1.2872  -0.4402 
## Trait.typeSurvivorship:Gradient.categoryResource     0.9904  -0.3994   0.4044 
##  
## Trait.typeFecundity:Gradient.categoryEnvironment 
## Trait.typeSurvivorship:Gradient.categoryEnvironment  *** 
## Trait.typeFecundity:Gradient.categoryPollution 
## Trait.typeSurvivorship:Gradient.categoryPollution      * 
## Trait.typeFecundity:Gradient.categoryResource        *** 
## Trait.typeSurvivorship:Gradient.categoryResource 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest.default(x= Q1$beta, sei =  Q1$se, ci.lb =  Q1$ci.lb, ci.ub =  Q1$ci.ub, 
                annotate=TRUE, showweights=T, header=F,
                slab = c("Environment fecundity", "Environment survivorship", "Pollution  fecundity", "Pollution  survivorship", "Resource  fecundity", "Resource  survivorship"))

Make Table 2 with results from Q1

Table2 <- data.frame(Stressor_type = c("Endogenous environment", "Endogenous environment", "Chemical pollution", "Chemical pollution", "Resource limitation","Resource limitation"),
                     Response_trait = c(rep(c("Fecundity", "Survivorship"), 3)),
                     Overall_mean = round(Q1$beta,3),
                     Lower_95 = round(Q1$ci.lb,3),
                     Upper_95 = round(Q1$ci.ub,3),
                     P_value = round(Q1$pval,3),
                     N_effects = c(length(dat_1$ID[dat_1$Gradient.category == "Environment" & dat_1$Trait.type == "Fecundity"]),
                                   length(dat_1$ID[dat_1$Gradient.category == "Environment" & dat_1$Trait.type == "Survivorship"]),
                                   length(dat_1$ID[dat_1$Gradient.category == "Pollution" & dat_1$Trait.type == "Fecundity"]),
                                   length(dat_1$ID[dat_1$Gradient.category == "Pollution" & dat_1$Trait.type == "Survivorship"]),
                                   length(dat_1$ID[dat_1$Gradient.category == "Resource" & dat_1$Trait.type == "Fecundity"]),
                                   length(dat_1$ID[dat_1$Gradient.category == "Resource" & dat_1$Trait.type == "Survivorship"])),
                     N_experiments = c(length(unique(dat_1$Experiment[dat_1$Gradient.category == "Environment" & dat_1$Trait.type == "Fecundity"])),
                                   length(unique(dat_1$Experiment[dat_1$Gradient.category == "Environment" & dat_1$Trait.type == "Survivorship"])),
                                   length(unique(dat_1$Experiment[dat_1$Gradient.category == "Pollution" & dat_1$Trait.type == "Fecundity"])),
                                   length(unique(dat_1$Experiment[dat_1$Gradient.category == "Pollution" & dat_1$Trait.type == "Survivorship"])),
                                   length(unique(dat_1$Experiment[dat_1$Gradient.category == "Resource" & dat_1$Trait.type == "Fecundity"])),
                                   length(unique(dat_1$Experiment[dat_1$Gradient.category == "Resource" & dat_1$Trait.type == "Survivorship"]))
),
N_host_taxa = c(length(unique(dat_1$Host[dat_1$Gradient.category == "Environment" & dat_1$Trait.type == "Fecundity"])),
                                   length(unique(dat_1$Host[dat_1$Gradient.category == "Environment" & dat_1$Trait.type == "Survivorship"])),
                                   length(unique(dat_1$Host[dat_1$Gradient.category == "Pollution" & dat_1$Trait.type == "Fecundity"])),
                                   length(unique(dat_1$Host[dat_1$Gradient.category == "Pollution" & dat_1$Trait.type == "Survivorship"])),
                                   length(unique(dat_1$Host[dat_1$Gradient.category == "Resource" & dat_1$Trait.type == "Fecundity"])),
                                   length(unique(dat_1$Host[dat_1$Gradient.category == "Resource" & dat_1$Trait.type == "Survivorship"]))
))

rownames(Table2) <- NULL

Table2 %>%
  kbl() %>%
  kable_material(c("striped", "hover"), full_width = F)
Stressor_type Response_trait Overall_mean Lower_95 Upper_95 P_value N_effects N_experiments N_host_taxa
Endogenous environment Fecundity -0.326 -0.802 0.150 0.180 14 3 3
Endogenous environment Survivorship -0.516 -0.715 -0.317 0.000 189 45 34
Chemical pollution Fecundity -0.141 -0.759 0.477 0.656 14 4 2
Chemical pollution Survivorship -0.351 -0.640 -0.062 0.017 76 23 15
Resource limitation Fecundity -0.864 -1.287 -0.440 0.000 59 6 4
Resource limitation Survivorship 0.002 -0.399 0.404 0.990 32 6 5

Make Figure 2 with the results of the best model for Q1

pal <- RColorBrewer::brewer.pal(11,"Spectral")

Q1_OP <-
  orchard_plot(
    Q1,
    xlab = "Standardised mean difference (Hedge's g)",
    transfm = "none",
    angle = 0,
    alpha = 0.7
  ) +
  scale_fill_manual(values = rep(pal[c(8,11,5)], each = 2)) + scale_colour_manual(values = rep(pal[c(8,11,5)], each = 2)) +
  theme(legend.position = c(0.25, 0.01), axis.title = element_text(size = 15), axis.text.y = element_text(size = 15)) +
  scale_y_discrete(
    labels = c("EE:fecundity",
               "EE:survivorship",
               "CP:fecundity",
               "CP:survivorship",
               "RL:fecundity",
               "RL:survivorship"))
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
Q1_OP

RVE

We check if our results are qualitatively similar if we use the robust variance estimator instead of the var-covar matrix for effects estimated in the same experiment.

Q1.R <-
  rma.mv(
    g ~ Trait.type:Gradient.category -1 ,
    V = dat_1$var.g,
    random = list( ~ 1 | Experiment, ~ 1 | ID),
    data = dat_1,
    method = "REML"
  )

Q1.R_est <- clubSandwich::coef_test(Q1.R, vcov = "CR2")
## Registered S3 method overwritten by 'clubSandwich':
##   method    from    
##   bread.mlm sandwich
forest.default(x= Q1.R_est$beta, sei =  Q1.R_est$SE, pch = 19,
                annotate=TRUE, showweights=F, header=c("Parameter", "Estimate [95% CI]"), 
                slab = c("EE:fecundity", "EE:survivorship", "CP:fecundity", "CP:survivorship", "RL:fecundity", "RL:survivorship"), 
               xlab = "Standardised mean difference (Hedge's g)")

Question 2 - Do fitness and infectivity traits have different sensitivity to stressors?

Model selection

All stressors

full_model.2 <- rma.mv(
    g,
    V = varcovmat_2_PD$mat,
    mods = ~ Trait.type*Gradient.category, 
    random = list( ~ 1 | Experiment,  ~ 1 | ID),
    data = dat_2,
    method = "ML"
  )

# Define if using Trait.type or Trait.category
model_selection.2 <- dredge(full_model.2, trace = 2) 
## Fixed term is "(Intercept)"
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=========                                                             |  12%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |==========================                                            |  38%
  |                                                                            
  |=============================================================         |  88%
subset(model_selection.2, delta <= 10, recalc.weights = FALSE)
## Global model call: rma.mv(yi = g, V = varcovmat_2_PD$mat, mods = ~Trait.type * Gradient.category, 
##     random = list(~1 | Experiment, ~1 | ID), data = dat_2, method = "ML")
## ---
## Model selection table 
##   (Int) Grd.ctg Trt.typ Grd.ctg:Trt.typ df    logLik   AICc delta weight
## 8     +       +       +               + 14 -1042.743 2114.1     0      1
## Models ranked by AICc(x)
# relative importance values for the predictors can be obtained with:
sw(model_selection.2)
##                      Trait.type Gradient.category Gradient.category:Trait.type
## Sum of weights:      1          1                 1                           
## N containing models: 3          3                 1

Fitness versus infectivity traits

full_model.2.rvf <- rma.mv(
    g,
    V = varcovmat_2_PD$mat,
    mods = ~ Trait.category*Gradient.category, 
    random = list( ~ 1 | Experiment,  ~ 1 | ID),
    data = dat_2,
    method = "ML"
  )
# Define if using Trait.type or Trait.category
model_selection.2.rvf <- dredge(full_model.2.rvf, trace = 2) 
## Fixed term is "(Intercept)"
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=========                                                             |  12%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |==========================                                            |  38%
  |                                                                            
  |=============================================================         |  88%
subset(model_selection.2.rvf, delta <= 10, recalc.weights = FALSE)
## Global model call: rma.mv(yi = g, V = varcovmat_2_PD$mat, mods = ~Trait.category * 
##     Gradient.category, random = list(~1 | Experiment, ~1 | ID), 
##     data = dat_2, method = "ML")
## ---
## Model selection table 
##   (Int) Grd.ctg Trt.ctg Grd.ctg:Trt.ctg df    logLik   AICc delta weight
## 8     +       +       +               +  8 -1053.678 2123.6     0      1
## Models ranked by AICc(x)
# relative importance values for the predictors can be obtained with:
sw(model_selection.2.rvf)
##                      Trait.category Gradient.category
## Sum of weights:      1              1                
## N containing models: 3              3                
##                      Gradient.category:Trait.category
## Sum of weights:      1                               
## N containing models: 1

Best model: Trait type x Gradient category interaction (lowest AICc)

Q2 <-
  rma.mv(
    g ~ Trait.type:Gradient.category -1,
    V = varcovmat_2_PD$mat,
    random = list(~ 1 | Experiment,  ~ 1 | ID),
    data = dat_2,
    method = "REML",
)

summary(Q2)
## 
## Multivariate Meta-Analysis Model (k = 686; method: REML)
## 
##     logLik    Deviance         AIC         BIC        AICc  ​ 
## -1031.4252   2062.8504   2090.8504   2154.0356   2091.4877   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed      factor 
## sigma^2.1  0.1697  0.4120    113     no  Experiment 
## sigma^2.2  0.4486  0.6698    686     no          ID 
## 
## Test for Residual Heterogeneity:
## QE(df = 674) = 1423699.2554, p-val < .0001
## 
## Test of Moderators (coefficients 1:12):
## QM(df = 12) = 103.7416, p-val < .0001
## 
## Model Results:
## 
##                                                      estimate      se     zval​ 
## Trait.typePrevalence:Gradient.categoryEnvironment      0.1010  0.1293   0.7813 
## Trait.typeIntensity:Gradient.categoryEnvironment       0.3978  0.0937   4.2434 
## Trait.typeFecundity:Gradient.categoryEnvironment       0.0318  0.2802   0.1134 
## Trait.typeSurvivorship:Gradient.categoryEnvironment   -0.5289  0.0901  -5.8734 
## Trait.typePrevalence:Gradient.categoryPollution       -0.5227  0.2520  -2.0741 
## Trait.typeIntensity:Gradient.categoryPollution        -0.0192  0.1266  -0.1518 
## Trait.typeFecundity:Gradient.categoryPollution         0.0887  0.2668   0.3324 
## Trait.typeSurvivorship:Gradient.categoryPollution     -0.2627  0.1318  -1.9933 
## Trait.typePrevalence:Gradient.categoryResource        -0.1503  0.2465  -0.6096 
## Trait.typeIntensity:Gradient.categoryResource         -0.4346  0.1826  -2.3803 
## Trait.typeFecundity:Gradient.categoryResource         -0.7416  0.2021  -3.6699 
## Trait.typeSurvivorship:Gradient.categoryResource       0.0049  0.1957   0.0253 
##                                                        pval    ci.lb    ci.ub 
## Trait.typePrevalence:Gradient.categoryEnvironment    0.4347  -0.1524   0.3545 
## Trait.typeIntensity:Gradient.categoryEnvironment     <.0001   0.2141   0.5815 
## Trait.typeFecundity:Gradient.categoryEnvironment     0.9097  -0.5173   0.5809 
## Trait.typeSurvivorship:Gradient.categoryEnvironment  <.0001  -0.7054  -0.3524 
## Trait.typePrevalence:Gradient.categoryPollution      0.0381  -1.0166  -0.0288 
## Trait.typeIntensity:Gradient.categoryPollution       0.8793  -0.2673   0.2289 
## Trait.typeFecundity:Gradient.categoryPollution       0.7396  -0.4343   0.6116 
## Trait.typeSurvivorship:Gradient.categoryPollution    0.0462  -0.5210  -0.0044 
## Trait.typePrevalence:Gradient.categoryResource       0.5421  -0.6335   0.3329 
## Trait.typeIntensity:Gradient.categoryResource        0.0173  -0.7924  -0.0767 
## Trait.typeFecundity:Gradient.categoryResource        0.0002  -1.1376  -0.3455 
## Trait.typeSurvivorship:Gradient.categoryResource     0.9799  -0.3786   0.3885 
##  
## Trait.typePrevalence:Gradient.categoryEnvironment 
## Trait.typeIntensity:Gradient.categoryEnvironment     *** 
## Trait.typeFecundity:Gradient.categoryEnvironment 
## Trait.typeSurvivorship:Gradient.categoryEnvironment  *** 
## Trait.typePrevalence:Gradient.categoryPollution        * 
## Trait.typeIntensity:Gradient.categoryPollution 
## Trait.typeFecundity:Gradient.categoryPollution 
## Trait.typeSurvivorship:Gradient.categoryPollution      * 
## Trait.typePrevalence:Gradient.categoryResource 
## Trait.typeIntensity:Gradient.categoryResource          * 
## Trait.typeFecundity:Gradient.categoryResource        *** 
## Trait.typeSurvivorship:Gradient.categoryResource 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest.default(x= Q2$beta, sei =  Q2$se, ci.lb =  Q2$ci.lb, ci.ub =  Q2$ci.ub, 
                annotate=TRUE, showweights=T, header=F,
                slab = c("EE:prevalence", "EE:intensity","EE:fecundity", "EE:survivorship",
               "CP:prevalence", "CP:intensity","CP:fecundity", "CP:survivorship",
               "RL:prevalence", "RL:intensity","RL:fecundity","RL:survivorship"))

par_short <- c("EE:prevalence", "EE:intensity","EE:fecundity", "EE:survivorship",
               "CP:prevalence", "CP:intensity","CP:fecundity", "CP:survivorship",
               "RL:prevalence", "RL:intensity","RL:fecundity","RL:survivorship")

pal <- RColorBrewer::brewer.pal(11,"Spectral")

Q2_OP <-
  orchard_plot(
    Q2,
    xlab = "Standardised mean difference (Hedge's g)",
    transfm = "none",
    angle = 0,
    alpha = 0.7
  ) +
  theme(legend.position = c(0.3, 0.01), axis.title = element_text(size = 15), axis.text.y = element_text(size = 15)) +
  scale_fill_manual(values = rep(pal[c(7,8,10,11,5,4)], each = 2)) +
  scale_colour_manual(values = rep(pal[c(7,8,10,11,5,4)], each = 2)) +
  scale_y_discrete(labels = par_short)
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
Q2_OP

Contrasts

Flipped values for fitness and infectivity Hedge’s g

dat_2$flipped <- dat_2$g
flip <- grep(pattern = "Intensity|Prevalence", x = dat_2$Trait.type)
for(i in flip){
  dat_2$flipped[i] <- -dat_2$flipped[i]
}
Re-run model with flipped values
Q2.flip <-
  rma.mv(
    flipped ~ Trait.type:Gradient.category -1,
    V = varcovmat_2_PD$mat,
    random = list(~ 1 | Experiment,  ~ 1 | ID),
    data = dat_2,
    method = "REML",
)

summary(Q2.flip)
## 
## Multivariate Meta-Analysis Model (k = 686; method: REML)
## 
##    logLik   Deviance        AIC        BIC       AICc  ​ 
## -979.7187  1959.4374  1987.4374  2050.6227  1988.0748   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed      factor 
## sigma^2.1  0.4654  0.6822    113     no  Experiment 
## sigma^2.2  0.3018  0.5493    686     no          ID 
## 
## Test for Residual Heterogeneity:
## QE(df = 674) = 1423699.2554, p-val < .0001
## 
## Test of Moderators (coefficients 1:12):
## QM(df = 12) = 67.2068, p-val < .0001
## 
## Model Results:
## 
##                                                      estimate      se     zval​ 
## Trait.typePrevalence:Gradient.categoryEnvironment     -0.3614  0.1363  -2.6525 
## Trait.typeIntensity:Gradient.categoryEnvironment      -0.4498  0.1086  -4.1415 
## Trait.typeFecundity:Gradient.categoryEnvironment      -0.4919  0.2574  -1.9108 
## Trait.typeSurvivorship:Gradient.categoryEnvironment   -0.6301  0.1047  -6.0168 
## Trait.typePrevalence:Gradient.categoryPollution        0.3523  0.2609   1.3503 
## Trait.typeIntensity:Gradient.categoryPollution         0.0600  0.1502   0.3995 
## Trait.typeFecundity:Gradient.categoryPollution        -0.2026  0.2611  -0.7759 
## Trait.typeSurvivorship:Gradient.categoryPollution     -0.2314  0.1525  -1.5172 
## Trait.typePrevalence:Gradient.categoryResource         0.1854  0.2365   0.7840 
## Trait.typeIntensity:Gradient.categoryResource          0.2959  0.1908   1.5506 
## Trait.typeFecundity:Gradient.categoryResource         -0.5041  0.2108  -2.3912 
## Trait.typeSurvivorship:Gradient.categoryResource      -0.2563  0.1940  -1.3214 
##                                                        pval    ci.lb    ci.ub 
## Trait.typePrevalence:Gradient.categoryEnvironment    0.0080  -0.6285  -0.0944 
## Trait.typeIntensity:Gradient.categoryEnvironment     <.0001  -0.6626  -0.2369 
## Trait.typeFecundity:Gradient.categoryEnvironment     0.0560  -0.9964   0.0127 
## Trait.typeSurvivorship:Gradient.categoryEnvironment  <.0001  -0.8354  -0.4248 
## Trait.typePrevalence:Gradient.categoryPollution      0.1769  -0.1591   0.8637 
## Trait.typeIntensity:Gradient.categoryPollution       0.6896  -0.2344   0.3543 
## Trait.typeFecundity:Gradient.categoryPollution       0.4378  -0.7145   0.3092 
## Trait.typeSurvivorship:Gradient.categoryPollution    0.1292  -0.5303   0.0675 
## Trait.typePrevalence:Gradient.categoryResource       0.4331  -0.2782   0.6490 
## Trait.typeIntensity:Gradient.categoryResource        0.1210  -0.0781   0.6699 
## Trait.typeFecundity:Gradient.categoryResource        0.0168  -0.9173  -0.0909 
## Trait.typeSurvivorship:Gradient.categoryResource     0.1864  -0.6365   0.1239 
##  
## Trait.typePrevalence:Gradient.categoryEnvironment     ** 
## Trait.typeIntensity:Gradient.categoryEnvironment     *** 
## Trait.typeFecundity:Gradient.categoryEnvironment       . 
## Trait.typeSurvivorship:Gradient.categoryEnvironment  *** 
## Trait.typePrevalence:Gradient.categoryPollution 
## Trait.typeIntensity:Gradient.categoryPollution 
## Trait.typeFecundity:Gradient.categoryPollution 
## Trait.typeSurvivorship:Gradient.categoryPollution 
## Trait.typePrevalence:Gradient.categoryResource 
## Trait.typeIntensity:Gradient.categoryResource 
## Trait.typeFecundity:Gradient.categoryResource          * 
## Trait.typeSurvivorship:Gradient.categoryResource 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest.default(x= Q2.flip$beta, sei =  Q2.flip$se, ci.lb =  Q2.flip$ci.lb, ci.ub =  Q2.flip$ci.ub, 
                annotate=TRUE, showweights=T, header=F,
                slab = c("EE:prevalence", "EE:intensity","EE:fecundity", "EE:survivorship",
               "CP:prevalence", "CP:intensity","CP:fecundity", "CP:survivorship",
               "RL:prevalence", "RL:intensity","RL:fecundity","RL:survivorship"))

Wald-type chi-square tests contrasting fitness vs infectivity effects on infected hosts with flipped values.
# Environment
# Prevalence vs Fecundity
anova(Q2.flip, X=c(-1,0,1,0,0,0,0,0,0,0,0,0))
## 
## Hypothesis:                                                                                                             
## 1: -Trait.typePrevalence:Gradient.categoryEnvironment + Trait.typeFecundity:Gradient.categoryEnvironment = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.1304 0.2672 -0.4882 0.6254 
## 
## Test of Hypothesis:
## QM(df = 1) = 0.2383, p-val = 0.6254
# Prevalence vs Survivorship
anova(Q2.flip, X=c(-1,0,0,1,0,0,0,0,0,0,0,0))
## 
## Hypothesis:                                                                                                                
## 1: -Trait.typePrevalence:Gradient.categoryEnvironment + Trait.typeSurvivorship:Gradient.categoryEnvironment = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.2687 0.1231 -2.1823 0.0291 
## 
## Test of Hypothesis:
## QM(df = 1) = 4.7626, p-val = 0.0291
# Intensity vs Fecundity
anova(Q2.flip, X=c(0,-1,1,0,0,0,0,0,0,0,0,0))
## 
## Hypothesis:                                                                                                            
## 1: -Trait.typeIntensity:Gradient.categoryEnvironment + Trait.typeFecundity:Gradient.categoryEnvironment = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.0421 0.2567 -0.1641 0.8697 
## 
## Test of Hypothesis:
## QM(df = 1) = 0.0269, p-val = 0.8697
# Intensity vs Survivorship
anova(Q2.flip, X=c(0,-1,0,1,0,0,0,0,0,0,0,0))
## 
## Hypothesis:                                                                                                               
## 1: -Trait.typeIntensity:Gradient.categoryEnvironment + Trait.typeSurvivorship:Gradient.categoryEnvironment = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.1803 0.0986 -1.8294 0.0673 
## 
## Test of Hypothesis:
## QM(df = 1) = 3.3468, p-val = 0.0673
# Pollution
# Prevalence vs Fecundity
anova(Q2.flip, X=c(0,0,0,0,-1,0,1,0,0,0,0,0))
## 
## Hypothesis:                                                                                                         
## 1: -Trait.typePrevalence:Gradient.categoryPollution + Trait.typeFecundity:Gradient.categoryPollution = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.5549 0.3376 -1.6437 0.1002 
## 
## Test of Hypothesis:
## QM(df = 1) = 2.7019, p-val = 0.1002
# Prevalence vs Survivorship
anova(Q2.flip, X=c(0,0,0,0,-1,0,0,1,0,0,0,0))
## 
## Hypothesis:                                                                                                            
## 1: -Trait.typePrevalence:Gradient.categoryPollution + Trait.typeSurvivorship:Gradient.categoryPollution = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.5837 0.2564 -2.2764 0.0228 
## 
## Test of Hypothesis:
## QM(df = 1) = 5.1819, p-val = 0.0228
# Intensity vs Fecundity
anova(Q2.flip, X=c(0,0,0,0,0,-1,1,0,0,0,0,0))
## 
## Hypothesis:                                                                                                        
## 1: -Trait.typeIntensity:Gradient.categoryPollution + Trait.typeFecundity:Gradient.categoryPollution = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.2626 0.2459 -1.0678 0.2856 
## 
## Test of Hypothesis:
## QM(df = 1) = 1.1402, p-val = 0.2856
# Intensity vs Survivorship
anova(Q2.flip, X=c(0,0,0,0,0,-1,0,1,0,0,0,0))
## 
## Hypothesis:                                                                                                           
## 1: -Trait.typeIntensity:Gradient.categoryPollution + Trait.typeSurvivorship:Gradient.categoryPollution = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.2914 0.1328 -2.1938 0.0283 
## 
## Test of Hypothesis:
## QM(df = 1) = 4.8126, p-val = 0.0283
# Resources
# Prevalence vs Fecundity
anova(Q2.flip, X=c(0,0,0,0,0,0,0,0,-1,0,1,0))
## 
## Hypothesis:                                                                                                       
## 1: -Trait.typePrevalence:Gradient.categoryResource + Trait.typeFecundity:Gradient.categoryResource = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.6895 0.2555 -2.6991 0.0070 
## 
## Test of Hypothesis:
## QM(df = 1) = 7.2852, p-val = 0.0070
# Prevalence vs Survivorship
anova(Q2.flip, X=c(0,0,0,0,0,0,0,0,-1,0,0,1))
## 
## Hypothesis:                                                                                                          
## 1: -Trait.typePrevalence:Gradient.categoryResource + Trait.typeSurvivorship:Gradient.categoryResource = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.4417 0.2596 -1.7013 0.0889 
## 
## Test of Hypothesis:
## QM(df = 1) = 2.8943, p-val = 0.0889
# Intensity vs Fecundity
anova(Q2.flip, X=c(0,0,0,0,0,0,0,0,0,-1,1,0))
## 
## Hypothesis:                                                                                                      
## 1: -Trait.typeIntensity:Gradient.categoryResource + Trait.typeFecundity:Gradient.categoryResource = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.8000 0.2025 -3.9509 <.0001 
## 
## Test of Hypothesis:
## QM(df = 1) = 15.6098, p-val < .0001
# Intensity vs Survivorship
anova(Q2.flip, X=c(0,0,0,0,0,0,0,0,0,-1,0,1))
## 
## Hypothesis:                                                                                                         
## 1: -Trait.typeIntensity:Gradient.categoryResource + Trait.typeSurvivorship:Gradient.categoryResource = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.5522 0.2128 -2.5950 0.0095 
## 
## Test of Hypothesis:
## QM(df = 1) = 6.7338, p-val = 0.0095

Absolute values for fitness and infectivity Hedge’s g

#absolute <- grep(pattern = "Fecundity|Survivorship|Intensity|Prevalence", x = dat_2$Trait.type)
#for(i in absolute){
#  dat_2$g[i] <- abs(dat_2$g[i])
#}

dat_2$absolute <- abs(dat_2$g)
Re-run model with absolute values
Q2.abs <-
  rma.mv(
    absolute ~ Trait.type:Gradient.category -1,
    V = varcovmat_2_PD$mat,
    random = list(~ 1 | Experiment,  ~ 1 | ID),
    data = dat_2,
    method = "REML",
)

summary(Q2.abs)
## 
## Multivariate Meta-Analysis Model (k = 686; method: REML)
## 
##    logLik   Deviance        AIC        BIC       AICc  ​ 
## -781.8300  1563.6600  1591.6600  1654.8452  1592.2973   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed      factor 
## sigma^2.1  0.1852  0.4304    113     no  Experiment 
## sigma^2.2  0.1073  0.3275    686     no          ID 
## 
## Test for Residual Heterogeneity:
## QE(df = 674) = 1427007.9678, p-val < .0001
## 
## Test of Moderators (coefficients 1:12):
## QM(df = 12) = 269.2974, p-val < .0001
## 
## Model Results:
## 
##                                                      estimate      se     zval​ 
## Trait.typePrevalence:Gradient.categoryEnvironment      0.8172  0.0997   8.1933 
## Trait.typeIntensity:Gradient.categoryEnvironment       0.9734  0.0771  12.6211 
## Trait.typeFecundity:Gradient.categoryEnvironment       0.5422  0.1850   2.9299 
## Trait.typeSurvivorship:Gradient.categoryEnvironment    0.7927  0.0738  10.7403 
## Trait.typePrevalence:Gradient.categoryPollution        0.7153  0.2012   3.5558 
## Trait.typeIntensity:Gradient.categoryPollution         0.6677  0.1039   6.4237 
## Trait.typeFecundity:Gradient.categoryPollution         0.8958  0.1966   4.5568 
## Trait.typeSurvivorship:Gradient.categoryPollution      0.6026  0.1060   5.6841 
## Trait.typePrevalence:Gradient.categoryResource         0.5347  0.1773   3.0152 
## Trait.typeIntensity:Gradient.categoryResource          0.8019  0.1376   5.8282 
## Trait.typeFecundity:Gradient.categoryResource          0.8126  0.1521   5.3440 
## Trait.typeSurvivorship:Gradient.categoryResource       0.8107  0.1406   5.7651 
##                                                        pval   ci.lb   ci.ub 
## Trait.typePrevalence:Gradient.categoryEnvironment    <.0001  0.6217  1.0127 
## Trait.typeIntensity:Gradient.categoryEnvironment     <.0001  0.8222  1.1245 
## Trait.typeFecundity:Gradient.categoryEnvironment     0.0034  0.1795  0.9048 
## Trait.typeSurvivorship:Gradient.categoryEnvironment  <.0001  0.6480  0.9374 
## Trait.typePrevalence:Gradient.categoryPollution      0.0004  0.3210  1.1096 
## Trait.typeIntensity:Gradient.categoryPollution       <.0001  0.4640  0.8714 
## Trait.typeFecundity:Gradient.categoryPollution       <.0001  0.5105  1.2811 
## Trait.typeSurvivorship:Gradient.categoryPollution    <.0001  0.3948  0.8103 
## Trait.typePrevalence:Gradient.categoryResource       0.0026  0.1871  0.8822 
## Trait.typeIntensity:Gradient.categoryResource        <.0001  0.5322  1.0715 
## Trait.typeFecundity:Gradient.categoryResource        <.0001  0.5145  1.1106 
## Trait.typeSurvivorship:Gradient.categoryResource     <.0001  0.5351  1.0863 
##  
## Trait.typePrevalence:Gradient.categoryEnvironment    *** 
## Trait.typeIntensity:Gradient.categoryEnvironment     *** 
## Trait.typeFecundity:Gradient.categoryEnvironment      ** 
## Trait.typeSurvivorship:Gradient.categoryEnvironment  *** 
## Trait.typePrevalence:Gradient.categoryPollution      *** 
## Trait.typeIntensity:Gradient.categoryPollution       *** 
## Trait.typeFecundity:Gradient.categoryPollution       *** 
## Trait.typeSurvivorship:Gradient.categoryPollution    *** 
## Trait.typePrevalence:Gradient.categoryResource        ** 
## Trait.typeIntensity:Gradient.categoryResource        *** 
## Trait.typeFecundity:Gradient.categoryResource        *** 
## Trait.typeSurvivorship:Gradient.categoryResource     *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest.default(x= Q2.abs$beta, sei =  Q2.abs$se, ci.lb =  Q2.abs$ci.lb, ci.ub =  Q2.abs$ci.ub, 
                annotate=TRUE, showweights=T, header=F,
                slab = c("EE:prevalence", "EE:intensity","EE:fecundity", "EE:survivorship",
               "CP:prevalence", "CP:intensity","CP:fecundity", "CP:survivorship",
               "RL:prevalence", "RL:intensity","RL:fecundity","RL:survivorship"))

Wald-type chi-square tests contrasting fitness vs infectivity effects on infected hosts.
# Environment
# Prevalence vs Fecundity
anova(Q2.abs, X=c(-1,0,1,0,0,0,0,0,0,0,0,0))
## 
## Hypothesis:                                                                                                             
## 1: -Trait.typePrevalence:Gradient.categoryEnvironment + Trait.typeFecundity:Gradient.categoryEnvironment = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.2751 0.1961 -1.4025 0.1608 
## 
## Test of Hypothesis:
## QM(df = 1) = 1.9670, p-val = 0.1608
# Prevalence vs Survivorship
anova(Q2.abs, X=c(-1,0,0,1,0,0,0,0,0,0,0,0))
## 
## Hypothesis:                                                                                                                
## 1: -Trait.typePrevalence:Gradient.categoryEnvironment + Trait.typeSurvivorship:Gradient.categoryEnvironment = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.0245 0.0946 -0.2592 0.7955 
## 
## Test of Hypothesis:
## QM(df = 1) = 0.0672, p-val = 0.7955
# Intensity vs Fecundity
anova(Q2.abs, X=c(0,-1,1,0,0,0,0,0,0,0,0,0))
## 
## Hypothesis:                                                                                                            
## 1: -Trait.typeIntensity:Gradient.categoryEnvironment + Trait.typeFecundity:Gradient.categoryEnvironment = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.4312 0.1874 -2.3008 0.0214 
## 
## Test of Hypothesis:
## QM(df = 1) = 5.2938, p-val = 0.0214
# Intensity vs Survivorship
anova(Q2.abs, X=c(0,-1,0,1,0,0,0,0,0,0,0,0))
## 
## Hypothesis:                                                                                                               
## 1: -Trait.typeIntensity:Gradient.categoryEnvironment + Trait.typeSurvivorship:Gradient.categoryEnvironment = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.1807 0.0778 -2.3232 0.0202 
## 
## Test of Hypothesis:
## QM(df = 1) = 5.3972, p-val = 0.0202
# Pollution
# Prevalence vs Fecundity
anova(Q2.abs, X=c(0,0,0,0,-1,0,1,0,0,0,0,0))
## 
## Hypothesis:                                                                                                         
## 1: -Trait.typePrevalence:Gradient.categoryPollution + Trait.typeFecundity:Gradient.categoryPollution = 0 
## 
## Results:
##    estimate     se   zval   pval 
## 1:   0.1805 0.2663 0.6778 0.4979 
## 
## Test of Hypothesis:
## QM(df = 1) = 0.4594, p-val = 0.4979
# Prevalence vs Survivorship
anova(Q2.abs, X=c(0,0,0,0,-1,0,0,1,0,0,0,0))
## 
## Hypothesis:                                                                                                            
## 1: -Trait.typePrevalence:Gradient.categoryPollution + Trait.typeSurvivorship:Gradient.categoryPollution = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.1127 0.2035 -0.5540 0.5795 
## 
## Test of Hypothesis:
## QM(df = 1) = 0.3070, p-val = 0.5795
# Intensity vs Fecundity
anova(Q2.abs, X=c(0,0,0,0,0,-1,1,0,0,0,0,0))
## 
## Hypothesis:                                                                                                        
## 1: -Trait.typeIntensity:Gradient.categoryPollution + Trait.typeFecundity:Gradient.categoryPollution = 0 
## 
## Results:
##    estimate     se   zval   pval 
## 1:   0.2281 0.1902 1.1993 0.2304 
## 
## Test of Hypothesis:
## QM(df = 1) = 1.4383, p-val = 0.2304
# Intensity vs Survivorship
anova(Q2.abs, X=c(0,0,0,0,0,-1,0,1,0,0,0,0))
## 
## Hypothesis:                                                                                                           
## 1: -Trait.typeIntensity:Gradient.categoryPollution + Trait.typeSurvivorship:Gradient.categoryPollution = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.0651 0.1023 -0.6363 0.5246 
## 
## Test of Hypothesis:
## QM(df = 1) = 0.4049, p-val = 0.5246
# Resources
# Prevalence vs Fecundity
anova(Q2.abs, X=c(0,0,0,0,0,0,0,0,-1,0,1,0))
## 
## Hypothesis:                                                                                                       
## 1: -Trait.typePrevalence:Gradient.categoryResource + Trait.typeFecundity:Gradient.categoryResource = 0 
## 
## Results:
##    estimate     se   zval   pval 
## 1:   0.2779 0.1922 1.4456 0.1483 
## 
## Test of Hypothesis:
## QM(df = 1) = 2.0898, p-val = 0.1483
# Prevalence vs Survivorship
anova(Q2.abs, X=c(0,0,0,0,0,0,0,0,-1,0,0,1))
## 
## Hypothesis:                                                                                                          
## 1: -Trait.typePrevalence:Gradient.categoryResource + Trait.typeSurvivorship:Gradient.categoryResource = 0 
## 
## Results:
##    estimate     se   zval   pval 
## 1:   0.2760 0.1994 1.3844 0.1662 
## 
## Test of Hypothesis:
## QM(df = 1) = 1.9165, p-val = 0.1662
# Intensity vs Fecundity
anova(Q2.abs, X=c(0,0,0,0,0,0,0,0,0,-1,1,0))
## 
## Hypothesis:                                                                                                      
## 1: -Trait.typeIntensity:Gradient.categoryResource + Trait.typeFecundity:Gradient.categoryResource = 0 
## 
## Results:
##    estimate     se   zval   pval 
## 1:   0.0107 0.1544 0.0692 0.9448 
## 
## Test of Hypothesis:
## QM(df = 1) = 0.0048, p-val = 0.9448
# Intensity vs Survivorship
anova(Q2.abs, X=c(0,0,0,0,0,0,0,0,0,-1,0,1))
## 
## Hypothesis:                                                                                                         
## 1: -Trait.typeIntensity:Gradient.categoryResource + Trait.typeSurvivorship:Gradient.categoryResource = 0 
## 
## Results:
##    estimate     se   zval   pval 
## 1:   0.0088 0.1604 0.0548 0.9563 
## 
## Test of Hypothesis:
## QM(df = 1) = 0.0030, p-val = 0.9563
Table3 <- data.frame(Stressor_type = rep(c("Endogenous environment",  "Chemical pollution","Resource limitation"), each = 4),
                     Response_trait = c(rep(c("Prevalence","Intensity","Fecundity", "Survivorship"), 3)),
                     Overall_mean = round(Q2$beta,3),
                     Lower_95 = round(Q2$ci.lb,3),
                     Upper_95 = round(Q2$ci.ub,3),
                     P_value = round(Q2$pval,3),
                     N_effects = c(length(dat_2$ID[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Prevalence"]),
                                   length(dat_2$ID[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Intensity"]),
                                   length(dat_2$ID[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Fecundity"]),
                                   length(dat_2$ID[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Survivorship"]),
                                   length(dat_2$ID[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Prevalence"]),
                                   length(dat_2$ID[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Intensity"]),
                                   length(dat_2$ID[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Fecundity"]),
                                   length(dat_2$ID[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Survivorship"]),
                                   length(dat_2$ID[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Prevalence"]),
                                   length(dat_2$ID[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Intensity"]),
                                   length(dat_2$ID[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Fecundity"]),
                                   length(dat_2$ID[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Survivorship"])),
                     N_experiments = c(length(unique(dat_2$Experiment[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Prevalence"])),
                                   length(unique(dat_2$Experiment[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Intensity"])),
                                   length(unique(dat_2$Experiment[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Fecundity"])),
                                   length(unique(dat_2$Experiment[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Survivorship"])),
                                   length(unique(dat_2$Experiment[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Prevalence"])),
                                   length(unique(dat_2$Experiment[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Intensity"])),
                                   length(unique(dat_2$Experiment[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Fecundity"])),
                                   length(unique(dat_2$Experiment[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Survivorship"])),
                                   length(unique(dat_2$Experiment[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Prevalence"])),
                                   length(unique(dat_2$Experiment[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Intensity"])),
                                   length(unique(dat_2$Experiment[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Fecundity"])),
                                   length(unique(dat_2$Experiment[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Survivorship"]))
),
N_host_taxa = c(length(unique(dat_2$Host[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Prevalence"])),
                length(unique(dat_2$Host[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Intensity"])),
                length(unique(dat_2$Host[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Fecundity"])),
                length(unique(dat_2$Host[dat_2$Gradient.category == "Environment" & dat_2$Trait.type == "Survivorship"])),
                length(unique(dat_2$Host[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Prevalence"])),
                length(unique(dat_2$Host[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Intensity"])),
                length(unique(dat_2$Host[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Fecundity"])),
                length(unique(dat_2$Host[dat_2$Gradient.category == "Pollution" & dat_2$Trait.type == "Survivorship"])),
                length(unique(dat_2$Host[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Prevalence"])),
                length(unique(dat_2$Host[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Intensity"])),
                length(unique(dat_2$Host[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Fecundity"])),
                length(unique(dat_2$Host[dat_2$Gradient.category == "Resource" & dat_2$Trait.type == "Survivorship"]))
))

rownames(Table3) <- NULL

Table3 %>%
  kbl() %>%
  kable_material(c("striped", "hover"), full_width = F)
Stressor_type Response_trait Overall_mean Lower_95 Upper_95 P_value N_effects N_experiments N_host_taxa
Endogenous environment Prevalence 0.101 -0.152 0.354 0.435 76 28 22
Endogenous environment Intensity 0.398 0.214 0.582 0.000 152 55 41
Endogenous environment Fecundity 0.032 -0.517 0.581 0.910 9 4 4
Endogenous environment Survivorship -0.529 -0.705 -0.352 0.000 162 66 47
Chemical pollution Prevalence -0.523 -1.017 -0.029 0.038 19 9 7
Chemical pollution Intensity -0.019 -0.267 0.229 0.879 74 30 20
Chemical pollution Fecundity 0.089 -0.434 0.612 0.740 16 7 4
Chemical pollution Survivorship -0.263 -0.521 -0.004 0.046 70 32 22
Resource limitation Prevalence -0.150 -0.633 0.333 0.542 15 5 4
Resource limitation Intensity -0.435 -0.792 -0.077 0.017 35 13 10
Resource limitation Fecundity -0.742 -1.138 -0.346 0.000 31 8 5
Resource limitation Survivorship 0.005 -0.379 0.388 0.980 27 10 9

RVE

We check if our results are qualitatively similar if we use the robust variance estimator instead of the var-covar matrix for effects estimated in the same experiment.

Q2.R <-
  rma.mv(
    g ~ Trait.type:Gradient.category -1,
    V = dat_2$var.g,
    random = list(~ 1 | Experiment,  ~ 1 | ID),
    data = dat_2,
    method = "REML",
)


Q2.R_est <- clubSandwich::coef_test(Q2.R, vcov = "CR2")

forest.default(x= Q2.R_est$beta, sei =  Q2.R_est$SE, pch = 19,
                annotate=TRUE, showweights=F, header=c("Parameter", "Estimate [95% CI]"), 
                slab = c("EE:prevalence", "EE:intensity", "EE:fecundity", "EE:survivorship", "CP:prevalence", "CP:intensity", "CP:fecundity", "CP:survivorship", "RL:prevalence", "RL:intensity", "RL:fecundity", "RL:survivorship"), 
               xlab = "Standardised mean difference (Hedge's g)")

Heterogeneity

Here I take the “best” model for each question to estimate between study heterogeneity. We could potentially do the same a different model or more models, depending on our goals.

Heterogeneity is important to know if the pooled effect size can be interpreted in a meaningful way. High heterogeneity is not necessarily a bad thing, but we would at least want to know if individual study effects usually go in the direction of the pooled effect.

I2 statistic

Is one of the most commonly reported measures of heterogeneity. It captures the “percentage of variability in the effect sizes that is not caused by sampling error”. That is, heterogeneity in true effect sizes. It tells us how much true heterogeneity there is relative to sampling error.

I2 formula

#' @title Parametric simulation 
#' @description Function for calculating I2 estimates using parametric simulations of model estimates taken from metafor. Note that the effectiveness of these simulations depends on the accuracy of model variance estimates.
#' @param estimate The estimate (i.e. variance) from a metafor model
#' @param sims The number of simulations 
#' @param n The sample size used in estimating the variance 
#' @author Daniel Noble - daniel.noble@anu.edu.au
#' @export
simMonteCarlo <- function(estimate, n, sims){
  tmp <- data.frame(num = base::rep(1:sims, each = n), 
                    y = stats::rnorm(n*sims, 0, base::sqrt(estimate)))
  Var <- tmp %>% dplyr::group_by(num) %>% dplyr::summarise(Mean_var = stats::var(y))
  return(as.numeric(Var$Mean_var))
}

## NOTE about PIPE: Run usethis::use_pipe() in the console. The package usethis will add what you need to import the pipe to your NAMESPACE and it will also drop warnings in checks

# Function for  rounding a data frame
round_df <- function(x, digits) {
  numeric_columns <- sapply(x, class) == 'numeric'
  x[numeric_columns] <-  round(x[numeric_columns], digits)
  x
}

# Function for estimating I2 
I2 <- function(model, v, ME = FALSE, sims = 1500, phylo = FALSE, obs = FALSE){
  
  if(class(model) != "rma.mv" && class(model) != "rma"){
    stop("The model object is not of class 'metafor'")
  }
  
  wi <- 1/v  #weight
  Vw <- sum((wi) * (length(wi) - 1))  / (((sum(wi)^2) - sum((wi)^2)))
  
  if("rma.mv" %in% class(model) | "rma" %in% class(model)){
    # Monte Carlo Simulations
    # From metafor extract the important statistics
    sigma2 <- matrix(model$sigma2, nrow = 1, ncol = length(model$sigma2))
    colnames(sigma2) <- model$s.names
    sigmaN <- model$s.nlevels
    
    if(obs == FALSE){
      stop("Please add the name of the observation-level random effect, obs. If models do not include this, re-run models including (~1|obs) in the random effect list")
    }
    
    #For each variance estimate use Monte Carlo simulation of data
    Sims <- data.frame(mapply(function(x,y) simMonteCarlo(x, y, sims = sims), x = sigma2, y = sigmaN))
    colnames(Sims) <- colnames(sigma2) 
    
    #Calculate total variance
    VT <- rowSums(cbind(Sims, Vw))
    Vt <- rowSums(Sims)  # remove Vw
    
    # For each variance component divide by the total variance. Note this needs to be fixed for phylo, but does deal with variable random effects.
    I2_re <- Sims / VT
    I2_total <- data.frame(Vt / VT)
    
    tmpMatrix <- data.frame(I2_re[, -match("ID", colnames(I2_re))], total = I2_total)
    names(tmpMatrix) = c(colnames(I2_re)[!colnames(I2_re) %in% 'ID'], 'total')
      
    CI <- lapply(tmpMatrix, function(x) stats::quantile(x, c(0.025, 0.975), na.rm = TRUE))
    I_CI <- as.data.frame(do.call(rbind, CI))
    colnames(I_CI) <- c("2.5% CI", "97.5% CI")
    I2_table <- cbind(I2_Est. = colMeans(tmpMatrix), I_CI )
    
    class(I2_table) <- c("metaAidR", "data.frame")
    
    return(round_df(I2_table, digits = 4))
  }
}

Question 1

First I calculate I2 using the formula above from Nakagawa and Santos (2012) for the model without the var-covar matrix

I2.Q1 <- I2(Q1.R, v = dat_1$var.g, obs = "ID")*100
## Warning in class(model) != "rma.mv" && class(model) != "rma": 'length(x) = 2 >
## 1' in coercion to 'logical(1)'
I2.Q1
##            I2_Est. 2.5% CI 97.5% CI
## Experiment   53.68   44.58    61.71
## total        93.86   92.71    94.86

Question 2

First I calculate I2 using the formula above from Nakagawa and Santos (2012) for the model without the var-covar matrix

I2.Q2 <- I2(Q2.R, v = dat_2$var.g, obs = "ID")*100
## Warning in class(model) != "rma.mv" && class(model) != "rma": 'length(x) = 2 >
## 1' in coercion to 'logical(1)'
I2.Q2
##            I2_Est. 2.5% CI 97.5% CI
## Experiment   26.22   20.97    31.31
## total        90.27   89.31    91.13

Publication bias

The most common way to assess publication bias is to test whether studies with large standard errors/low precision/small sample size tend to have larger effect sizes. Small studies have larger SEs, and because they are small, they are more likely to get published only if results are significant. In contrast, large studies, which have higher precision and have costed more money, are more likely to get published either way. For this reason small studies can cause publication bias.

Funnel plots

A very common way to assess publication bias is therefore to check whether studies with low precision generally have larger effects. The typical way to visualize this is through a funnel plot with all effects. Now, this is not entirely correct for us, as the effects are not independent. Nakagawa et al 2021 (Methods Ecol Evol) suggest using conditional residuals for funnel plots. These residual plots take some of the non-independence into account, but are not perfect, as they assume: 1)sampling SE(sei) does not covary with moderators in meta-regression, 2) sampling SE is the same as the SE of the residuals and 3) sampling variances are not correlated. This last one is most certainly incorrect.

Estimation of conditional residuals is not available to rma.mv objects (which can include a sampling variance-covariance matrix). The available alternative is to ignore this covariation in sampling variance for the assessment of publication bias. This is a caveat that should be acknowledgend in the paper, yet I think it is a useful visualization nonetheless. To do these funnel plots I re-fit the best model as rma.uni for each data set. While not ideal, this the best current alternative suggested by Nakagawa et al 2021 (Methods Ecol Evol)

Note I checked some previous studies that use multilevel regression models as we do here, to see what they do about publication bias. Sauer et al 2020 Ecology (101) acknowledged this problem and decided not to even show funnel plots: “Because of the complex non independence among effect sizes within a study (e.g., some studies had multiple effect sizes), we did not use funnel plots or rank correlation tests to assess publication bias (Lau et al. 2006, Civitello et al.2015).” Others (Yates et al. 2019 Environmental DNA (1); Salerno et al 2021 Global Change Biology (27); Brlik et al. 2020 JAE) acknowledge this problem, but still show funnel plots warning caution. They also conduct a version of Eggerts test somewhat in the same spirit as I do below.

Data set 1: fitness effects of environmental stress

# re fit as a mixed effects model with rma.uni
# any fitting method other than "FE" works for a random/mixed effects model
Q1.uni <-
  rma.uni(
    g ~ Trait.type:Gradient.category + Experiment + ID ,
    vi = dat_1$var.g,
    data = dat_1,
    method = "REML"
  )
## Warning: Redundant predictors dropped from the model.
# get conditional residuals (residuals for each effect)
Q1res <- rstandard.rma.uni(Q1.uni, type = "conditional")

# create colour palette as  gradient of effect size g
nlvl <- length(dat_1$g)
pal <- colorRampPalette(brewer.pal(9, "YlOrRd"))(nlvl)
pal <-pal[order(match(pal,dat_1$g))]

# plot funnel with conditional residuals
metafor::funnel(Q1res$resid, sqrt(dat_1$var.g), level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), col=pal, back = "white", refline=0, legend=TRUE, xlab = "Conditional residual", ylab = "Standard error")

We can see that as SE increases (e.g ~ 1) all studies report large negative fitness effects of environmental stress (large negative residuals). In other words, it seems that very small studies tend to find that environmental stress is more negative to fitness!

Data set 2: Fitness and infectivity effects of environmental stress

# re fit as a mixed effects model with rma.uni
# any fitting method other than "FE" works for a random/mixed effects model
Q2.uni <-
  rma.uni(
    g ~ Trait.type:Gradient.category + Experiment + ID,
    vi = dat_2$var.g,
    data = dat_2,
    method = "REML"
  )
## Warning: Redundant predictors dropped from the model.
# get conditional residuals (residuals for each effect)
Q2res <- rstandard.rma.uni(Q2.uni, type = "conditional")

# create colour palette as  gradient of effect size g
nlvl <- length(dat_2$g)
pal <- colorRampPalette(brewer.pal(9, "YlOrRd"))(nlvl)
pal <-pal[order(match(pal,dat_2$g))]

# plot funnel with all conditional residuals
metafor::funnel(Q2res$resid, sqrt(dat_2$var.g), level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), col=pal, back = "white", refline=0, legend=TRUE, xlab = "Conditional residual", ylab = "Standard error")

Because this data set includes both the infectivity and fitness effects it’s not surprising that small studies are distributed on both sides of the funnel. This does not mean absence of publication bias as negative effects for fitness and positive effects for (-) infectivity both would suggest that small studies tend to inflate the negative consequences of environmental stress. For the publication, I suggest we use -g for infectivity effects (currently g = infection intensity or prevalence) and/or we use funnel plots for only dataset 1.

Two step approach to test for publication bias

Following Nakagawa et al. (2021) we conduct a two-step approach to test for publication bias, that can be thought of as a modified Egger’s test for multilevel meta analysis. First we include the effect sizes’ standard error as the only moderator and including the same random effect structure as in our meta-analysis. If the slope of this moderator is significant, the estimated intercept provides a downwardly biased estimate of the overall meta-analytic effect (Doucouliagos 2012, 2014). In this case, we run a second meta regression with only the effect sizes’ variances as moderators. The intercept in this case may still be biased, but less so.

Data set 1: fitness effects of environmental stress

# get standard errors of effect sizes 
dat_1$sei <- sqrt(dat_1$var.g)

# meta regression with SE
d1s1 <- rma.mv(
  g ~ 1 + sei,
  V = varcovmat_1_PD$mat,
  random = list(~ 1 | Experiment, ~ 1 | ID),
  data = dat_1,
  method = "REML"
)

d1s1
## 
## Multivariate Meta-Analysis Model (k = 384; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed      factor 
## sigma^2.1  0.4132  0.6428     78     no  Experiment 
## sigma^2.2  0.2912  0.5396    384     no          ID 
## 
## Test for Residual Heterogeneity:
## QE(df = 382) = 247179.2283, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## QM(df = 1) = 56.3837, p-val < .0001
## 
## Model Results:
## 
##          estimate      se     zval    pval    ci.lb    ci.ub     ​ 
## intrcpt    0.3606  0.1412   2.5528  0.0107   0.0837   0.6374    * 
## sei       -1.5464  0.2059  -7.5089  <.0001  -1.9501  -1.1428  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The slope is negative and significant. Small studies report large and negative effect sizes!! Because the slope of sei is significant, variance is a better moderator according to (Doucouliagos 2012, 2014).

d1s2 <- rma.mv(
  g ~ 1 + dat_1$var.g,
  V = varcovmat_1_PD$mat,
  random = list(~ 1 | Experiment, ~ 1 | ID),
  data = dat_1,
  method = "REML"
)

d1s2 
## 
## Multivariate Meta-Analysis Model (k = 384; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed      factor 
## sigma^2.1  0.3473  0.5893     78     no  Experiment 
## sigma^2.2  0.2943  0.5425    384     no          ID 
## 
## Test for Residual Heterogeneity:
## QE(df = 382) = 422619.1025, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## QM(df = 1) = 77.9257, p-val < .0001
## 
## Model Results:
## 
##              estimate      se     zval    pval    ci.lb    ci.ub     ​ 
## intrcpt       -0.1613  0.0911  -1.7708  0.0766  -0.3398   0.0172    . 
## dat_1$var.g   -0.8711  0.0987  -8.8276  <.0001  -1.0645  -0.6777  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Again, precision effects are significant, suggesting publication bias. Thus, to “account” for potential effects of publication bias, we need to include variance as a moderator.

Including the variance-covariance error matrix within studies

Q1.B <-
  rma.mv(
    g ~ Trait.type:Gradient.category + var.g -1,
    V = varcovmat_1_PD$mat,
    random = list( ~ 1 | Experiment, ~ 1 | ID),
    data = dat_1,
    method = "REML"
  )

summary(Q1.B)
## 
## Multivariate Meta-Analysis Model (k = 384; method: REML)
## 
##    logLik   Deviance        AIC        BIC       AICc  ​ 
## -507.8909  1015.7818  1033.7818  1069.1720  1034.2723   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed      factor 
## sigma^2.1  0.3149  0.5612     78     no  Experiment 
## sigma^2.2  0.2826  0.5316    384     no          ID 
## 
## Test for Residual Heterogeneity:
## QE(df = 377) = 133298.1893, p-val < .0001
## 
## Test of Moderators (coefficients 1:7):
## QM(df = 7) = 125.9072, p-val < .0001
## 
## Model Results:
## 
##                                                      estimate      se     zval​ 
## var.g                                                 -0.8846  0.0980  -9.0217 
## Trait.typeFecundity:Gradient.categoryEnvironment      -0.0978  0.2332  -0.4196 
## Trait.typeSurvivorship:Gradient.categoryEnvironment   -0.1762  0.1109  -1.5890 
## Trait.typeFecundity:Gradient.categoryPollution         0.0589  0.3121   0.1886 
## Trait.typeSurvivorship:Gradient.categoryPollution     -0.1061  0.1548  -0.6853 
## Trait.typeFecundity:Gradient.categoryResource         -0.6939  0.2177  -3.1877 
## Trait.typeSurvivorship:Gradient.categoryResource       0.3048  0.2032   1.5002 
##                                                        pval    ci.lb    ci.ub 
## var.g                                                <.0001  -1.0767  -0.6924 
## Trait.typeFecundity:Gradient.categoryEnvironment     0.6747  -0.5548   0.3591 
## Trait.typeSurvivorship:Gradient.categoryEnvironment  0.1121  -0.3934   0.0411 
## Trait.typeFecundity:Gradient.categoryPollution       0.8504  -0.5529   0.6706 
## Trait.typeSurvivorship:Gradient.categoryPollution    0.4932  -0.4094   0.1973 
## Trait.typeFecundity:Gradient.categoryResource        0.0014  -1.1206  -0.2673 
## Trait.typeSurvivorship:Gradient.categoryResource     0.1336  -0.0934   0.7030 
##  
## var.g                                                *** 
## Trait.typeFecundity:Gradient.categoryEnvironment 
## Trait.typeSurvivorship:Gradient.categoryEnvironment 
## Trait.typeFecundity:Gradient.categoryPollution 
## Trait.typeSurvivorship:Gradient.categoryPollution 
## Trait.typeFecundity:Gradient.categoryResource         ** 
## Trait.typeSurvivorship:Gradient.categoryResource 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest.default(x= Q1.B$beta, sei =  Q1.B$se, pch = 19,
                annotate=TRUE, showweights=F, header=c("Parameter", "Estimate [95% CI]"), 
                slab = c("variance", "EE:fecundity", "EE:survivorship", "CP:fecundity", "CP:survivorship", "RL:fecundity", "RL:survivorship"), 
               xlab = "Standardised mean difference (Hedge's g)")

Data set 2: fitness and infectivity effects of environmental stress

# get standard errors of effect sizes 
dat_2$sei <- sqrt(dat_2$var.g)

# meta regression with SE
d2s1 <- rma.mv(
  g ~ 1 + sei,
  V = varcovmat_2_PD$mat,
  random = list(~ 1 | Experiment, ~ 1 | ID),
  data = dat_2,
  method = "REML"
)

As in data set 1, effects of small studies are dominated by negative effect sizes. Because the slope of sei is significant the variance is a better moderator for the overall effect.

# meta regression with SE
d2s2 <- rma.mv(
  g ~ 1 + dat_2$var.g,
  V = varcovmat_2_PD$mat,
  random = list(~ 1 | Experiment, ~ 1 | ID),
  data = dat_2,
  method = "REML"
)

d2s2
## 
## Multivariate Meta-Analysis Model (k = 686; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed      factor 
## sigma^2.1  0.1407  0.3751    113     no  Experiment 
## sigma^2.2  0.5141  0.7170    686     no          ID 
## 
## Test for Residual Heterogeneity:
## QE(df = 684) = 3689465.0418, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## QM(df = 1) = 16.9025, p-val < .0001
## 
## Model Results:
## 
##              estimate      se     zval    pval    ci.lb    ci.ub     ​ 
## intrcpt       -0.0081  0.0597  -0.1363  0.8916  -0.1251   0.1088      
## dat_2$var.g   -0.3665  0.0891  -4.1113  <.0001  -0.5411  -0.1918  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Q2.B <-
  rma.mv(
    g ~ Trait.type:Gradient.category + var.g -1,
    V = varcovmat_2_PD$mat,
    random = list( ~ 1 | Experiment, ~ 1 | ID),
    data = dat_2,
    method = "REML"
  )

summary(Q2.B)
## 
## Multivariate Meta-Analysis Model (k = 686; method: REML)
## 
##     logLik    Deviance         AIC         BIC        AICc  ​ 
## -1020.1000   2040.2001   2070.2001   2137.8762   2070.9307   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed      factor 
## sigma^2.1  0.1691  0.4112    113     no  Experiment 
## sigma^2.2  0.4257  0.6525    686     no          ID 
## 
## Test for Residual Heterogeneity:
## QE(df = 673) = 827885.4316, p-val < .0001
## 
## Test of Moderators (coefficients 1:13):
## QM(df = 13) = 124.3757, p-val < .0001
## 
## Model Results:
## 
##                                                      estimate      se     zval​ 
## var.g                                                 -0.3826  0.0884  -4.3297 
## Trait.typePrevalence:Gradient.categoryEnvironment      0.2548  0.1326   1.9217 
## Trait.typeIntensity:Gradient.categoryEnvironment       0.5017  0.0958   5.2366 
## Trait.typeFecundity:Gradient.categoryEnvironment       0.1025  0.2754   0.3723 
## Trait.typeSurvivorship:Gradient.categoryEnvironment   -0.4050  0.0932  -4.3449 
## Trait.typePrevalence:Gradient.categoryPollution       -0.3950  0.2506  -1.5764 
## Trait.typeIntensity:Gradient.categoryPollution         0.0501  0.1261   0.3970 
## Trait.typeFecundity:Gradient.categoryPollution         0.2032  0.2645   0.7685 
## Trait.typeSurvivorship:Gradient.categoryPollution     -0.1509  0.1326  -1.1378 
## Trait.typePrevalence:Gradient.categoryResource        -0.0399  0.2440  -0.1635 
## Trait.typeIntensity:Gradient.categoryResource         -0.3309  0.1816  -1.8224 
## Trait.typeFecundity:Gradient.categoryResource         -0.6764  0.1999  -3.3847 
## Trait.typeSurvivorship:Gradient.categoryResource       0.1080  0.1942   0.5562 
##                                                        pval    ci.lb    ci.ub 
## var.g                                                <.0001  -0.5557  -0.2094 
## Trait.typePrevalence:Gradient.categoryEnvironment    0.0546  -0.0051   0.5146 
## Trait.typeIntensity:Gradient.categoryEnvironment     <.0001   0.3139   0.6895 
## Trait.typeFecundity:Gradient.categoryEnvironment     0.7097  -0.4373   0.6424 
## Trait.typeSurvivorship:Gradient.categoryEnvironment  <.0001  -0.5877  -0.2223 
## Trait.typePrevalence:Gradient.categoryPollution      0.1149  -0.8862   0.0961 
## Trait.typeIntensity:Gradient.categoryPollution       0.6913  -0.1971   0.2973 
## Trait.typeFecundity:Gradient.categoryPollution       0.4422  -0.3151   0.7216 
## Trait.typeSurvivorship:Gradient.categoryPollution    0.2552  -0.4108   0.1090 
## Trait.typePrevalence:Gradient.categoryResource       0.8701  -0.5182   0.4384 
## Trait.typeIntensity:Gradient.categoryResource        0.0684  -0.6868   0.0250 
## Trait.typeFecundity:Gradient.categoryResource        0.0007  -1.0681  -0.2847 
## Trait.typeSurvivorship:Gradient.categoryResource     0.5781  -0.2726   0.4886 
##  
## var.g                                                *** 
## Trait.typePrevalence:Gradient.categoryEnvironment      . 
## Trait.typeIntensity:Gradient.categoryEnvironment     *** 
## Trait.typeFecundity:Gradient.categoryEnvironment 
## Trait.typeSurvivorship:Gradient.categoryEnvironment  *** 
## Trait.typePrevalence:Gradient.categoryPollution 
## Trait.typeIntensity:Gradient.categoryPollution 
## Trait.typeFecundity:Gradient.categoryPollution 
## Trait.typeSurvivorship:Gradient.categoryPollution 
## Trait.typePrevalence:Gradient.categoryResource 
## Trait.typeIntensity:Gradient.categoryResource          . 
## Trait.typeFecundity:Gradient.categoryResource        *** 
## Trait.typeSurvivorship:Gradient.categoryResource 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest.default(x= Q2.B$beta, sei =  Q2.B$se, pch = 19,
                annotate=TRUE, showweights=F, header=c("Parameter", "Estimate [95% CI]"), 
                slab = c("variance","EE:prevalence", "EE:intensity", "EE:fecundity", "EE:survivorship", "CP:prevalence", "CP:intensity", "CP:fecundity", "CP:survivorship", "RL:prevalence", "RL:intensity", "RL:fecundity", "RL:survivorship"), 
               xlab = "Standardised mean difference (Hedge's g)")

Contrasts

Re-run model with absolute values
Q2.B.abs <-
  rma.mv(
    absolute ~ Trait.type:Gradient.category + var.g -1,
    V = varcovmat_2_PD$mat,
    random = list(~ 1 | Experiment,  ~ 1 | ID),
    data = dat_2,
    method = "REML",
)

summary(Q2.B.abs)
## 
## Multivariate Meta-Analysis Model (k = 686; method: REML)
## 
##    logLik   Deviance        AIC        BIC       AICc  ​ 
## -671.1227  1342.2454  1372.2454  1439.9216  1372.9760   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed      factor 
## sigma^2.1  0.0890  0.2983    113     no  Experiment 
## sigma^2.2  0.0712  0.2669    686     no          ID 
## 
## Test for Residual Heterogeneity:
## QE(df = 673) = 822493.8021, p-val < .0001
## 
## Test of Moderators (coefficients 1:13):
## QM(df = 13) = 633.2958, p-val < .0001
## 
## Model Results:
## 
##                                                      estimate      se     zval​ 
## var.g                                                  1.0912  0.0709  15.4002 
## Trait.typePrevalence:Gradient.categoryEnvironment      0.4348  0.0888   4.8977 
## Trait.typeIntensity:Gradient.categoryEnvironment       0.6922  0.0661  10.4710 
## Trait.typeFecundity:Gradient.categoryEnvironment       0.3123  0.1641   1.9036 
## Trait.typeSurvivorship:Gradient.categoryEnvironment    0.5140  0.0625   8.2278 
## Trait.typePrevalence:Gradient.categoryPollution        0.3522  0.1796   1.9608 
## Trait.typeIntensity:Gradient.categoryPollution         0.4506  0.0850   5.2989 
## Trait.typeFecundity:Gradient.categoryPollution         0.5483  0.1760   3.1149 
## Trait.typeSurvivorship:Gradient.categoryPollution      0.3254  0.0879   3.7015 
## Trait.typePrevalence:Gradient.categoryResource         0.2295  0.1593   1.4409 
## Trait.typeIntensity:Gradient.categoryResource          0.4900  0.1185   4.1343 
## Trait.typeFecundity:Gradient.categoryResource          0.5740  0.1297   4.4262 
## Trait.typeSurvivorship:Gradient.categoryResource       0.5195  0.1231   4.2215 
##                                                        pval    ci.lb   ci.ub 
## var.g                                                <.0001   0.9524  1.2301 
## Trait.typePrevalence:Gradient.categoryEnvironment    <.0001   0.2608  0.6088 
## Trait.typeIntensity:Gradient.categoryEnvironment     <.0001   0.5626  0.8218 
## Trait.typeFecundity:Gradient.categoryEnvironment     0.0570  -0.0093  0.6338 
## Trait.typeSurvivorship:Gradient.categoryEnvironment  <.0001   0.3916  0.6365 
## Trait.typePrevalence:Gradient.categoryPollution      0.0499   0.0002  0.7041 
## Trait.typeIntensity:Gradient.categoryPollution       <.0001   0.2839  0.6172 
## Trait.typeFecundity:Gradient.categoryPollution       0.0018   0.2033  0.8933 
## Trait.typeSurvivorship:Gradient.categoryPollution    0.0002   0.1531  0.4977 
## Trait.typePrevalence:Gradient.categoryResource       0.1496  -0.0827  0.5417 
## Trait.typeIntensity:Gradient.categoryResource        <.0001   0.2577  0.7223 
## Trait.typeFecundity:Gradient.categoryResource        <.0001   0.3198  0.8281 
## Trait.typeSurvivorship:Gradient.categoryResource     <.0001   0.2783  0.7607 
##  
## var.g                                                *** 
## Trait.typePrevalence:Gradient.categoryEnvironment    *** 
## Trait.typeIntensity:Gradient.categoryEnvironment     *** 
## Trait.typeFecundity:Gradient.categoryEnvironment       . 
## Trait.typeSurvivorship:Gradient.categoryEnvironment  *** 
## Trait.typePrevalence:Gradient.categoryPollution        * 
## Trait.typeIntensity:Gradient.categoryPollution       *** 
## Trait.typeFecundity:Gradient.categoryPollution        ** 
## Trait.typeSurvivorship:Gradient.categoryPollution    *** 
## Trait.typePrevalence:Gradient.categoryResource 
## Trait.typeIntensity:Gradient.categoryResource        *** 
## Trait.typeFecundity:Gradient.categoryResource        *** 
## Trait.typeSurvivorship:Gradient.categoryResource     *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest.default(x= Q2.B.abs$beta, sei =  Q2.B.abs$se, pch = 19,
                annotate=TRUE, showweights=F, header=c("Parameter", "Estimate [95% CI]"), 
                slab = c("variance","EE:prevalence", "EE:intensity", "EE:fecundity", "EE:survivorship", "CP:prevalence", "CP:intensity", "CP:fecundity", "CP:survivorship", "RL:prevalence", "RL:intensity", "RL:fecundity", "RL:survivorship"), 
               xlab = "Standardised mean difference (Hedge's g)")

Wald-type chi-square tests contrasting fitness vs infectivity effects on infected hosts.

# Environment
# Prevalence vs Fecundity
anova(Q2.B.abs, X=c(0,-1,0,1,0,0,0,0,0,0,0,0,0))
## 
## Hypothesis:                                                                                                             
## 1: -Trait.typePrevalence:Gradient.categoryEnvironment + Trait.typeFecundity:Gradient.categoryEnvironment = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.1225 0.1768 -0.6930 0.4883 
## 
## Test of Hypothesis:
## QM(df = 1) = 0.4802, p-val = 0.4883
# Prevalence vs Survivorship
anova(Q2.B.abs, X=c(0,-1,0,0,1,0,0,0,0,0,0,0,0))
## 
## Hypothesis:                                                                                                                
## 1: -Trait.typePrevalence:Gradient.categoryEnvironment + Trait.typeSurvivorship:Gradient.categoryEnvironment = 0 
## 
## Results:
##    estimate     se   zval   pval 
## 1:   0.0792 0.0865 0.9156 0.3599 
## 
## Test of Hypothesis:
## QM(df = 1) = 0.8383, p-val = 0.3599
# Intensity vs Fecundity
anova(Q2.B.abs, X=c(0,0,-1,1,0,0,0,0,0,0,0,0,0))
## 
## Hypothesis:                                                                                                            
## 1: -Trait.typeIntensity:Gradient.categoryEnvironment + Trait.typeFecundity:Gradient.categoryEnvironment = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.3799 0.1684 -2.2559 0.0241 
## 
## Test of Hypothesis:
## QM(df = 1) = 5.0889, p-val = 0.0241
# Intensity vs Survivorship
anova(Q2.B.abs, X=c(0,0,-1,0,1,0,0,0,0,0,0,0,0))
## 
## Hypothesis:                                                                                                               
## 1: -Trait.typeIntensity:Gradient.categoryEnvironment + Trait.typeSurvivorship:Gradient.categoryEnvironment = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.1782 0.0716 -2.4885 0.0128 
## 
## Test of Hypothesis:
## QM(df = 1) = 6.1928, p-val = 0.0128
# Pollution
# Prevalence vs Fecundity
anova(Q2.B.abs, X=c(0,0,0,0,0,-1,0,1,0,0,0,0,0))
## 
## Hypothesis:                                                                                                         
## 1: -Trait.typePrevalence:Gradient.categoryPollution + Trait.typeFecundity:Gradient.categoryPollution = 0 
## 
## Results:
##    estimate     se   zval   pval 
## 1:   0.1962 0.2431 0.8070 0.4197 
## 
## Test of Hypothesis:
## QM(df = 1) = 0.6512, p-val = 0.4197
# Prevalence vs Survivorship
anova(Q2.B.abs, X=c(0,0,0,0,0,-1,0,0,1,0,0,0,0))
## 
## Hypothesis:                                                                                                            
## 1: -Trait.typePrevalence:Gradient.categoryPollution + Trait.typeSurvivorship:Gradient.categoryPollution = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.0267 0.1856 -0.1441 0.8854 
## 
## Test of Hypothesis:
## QM(df = 1) = 0.0208, p-val = 0.8854
# Intensity vs Fecundity
anova(Q2.B.abs, X=c(0,0,0,0,0,0,-1,1,0,0,0,0,0))
## 
## Hypothesis:                                                                                                        
## 1: -Trait.typeIntensity:Gradient.categoryPollution + Trait.typeFecundity:Gradient.categoryPollution = 0 
## 
## Results:
##    estimate     se   zval   pval 
## 1:   0.0977 0.1747 0.5593 0.5759 
## 
## Test of Hypothesis:
## QM(df = 1) = 0.3129, p-val = 0.5759
# Intensity vs Survivorship
anova(Q2.B.abs, X=c(0,0,0,0,0,0,-1,0,1,0,0,0,0))
## 
## Hypothesis:                                                                                                           
## 1: -Trait.typeIntensity:Gradient.categoryPollution + Trait.typeSurvivorship:Gradient.categoryPollution = 0 
## 
## Results:
##    estimate     se    zval   pval 
## 1:  -0.1252 0.0937 -1.3358 0.1816 
## 
## Test of Hypothesis:
## QM(df = 1) = 1.7843, p-val = 0.1816
# Resources
# Prevalence vs Fecundity
anova(Q2.B.abs, X=c(0,0,0,0,0,0,0,0,0,-1,0,1,0))
## 
## Hypothesis:                                                                                                       
## 1: -Trait.typePrevalence:Gradient.categoryResource + Trait.typeFecundity:Gradient.categoryResource = 0 
## 
## Results:
##    estimate     se   zval   pval 
## 1:   0.3445 0.1754 1.9640 0.0495 
## 
## Test of Hypothesis:
## QM(df = 1) = 3.8574, p-val = 0.0495
# Prevalence vs Survivorship
anova(Q2.B.abs, X=c(0,0,0,0,0,0,0,0,0,-1,0,0,1))
## 
## Hypothesis:                                                                                                          
## 1: -Trait.typePrevalence:Gradient.categoryResource + Trait.typeSurvivorship:Gradient.categoryResource = 0 
## 
## Results:
##    estimate     se   zval   pval 
## 1:   0.2900 0.1826 1.5885 0.1122 
## 
## Test of Hypothesis:
## QM(df = 1) = 2.5233, p-val = 0.1122
# Intensity vs Fecundity
anova(Q2.B.abs, X=c(0,0,0,0,0,0,0,0,0,0,-1,1,0))
## 
## Hypothesis:                                                                                                      
## 1: -Trait.typeIntensity:Gradient.categoryResource + Trait.typeFecundity:Gradient.categoryResource = 0 
## 
## Results:
##    estimate     se   zval   pval 
## 1:   0.0840 0.1405 0.5975 0.5502 
## 
## Test of Hypothesis:
## QM(df = 1) = 0.3570, p-val = 0.5502
# Intensity vs Survivorship
anova(Q2.B.abs, X=c(0,0,0,0,0,0,0,0,0,0,-1,0,1))
## 
## Hypothesis:                                                                                                         
## 1: -Trait.typeIntensity:Gradient.categoryResource + Trait.typeSurvivorship:Gradient.categoryResource = 0 
## 
## Results:
##    estimate     se   zval   pval 
## 1:   0.0295 0.1452 0.2031 0.8391 
## 
## Test of Hypothesis:
## QM(df = 1) = 0.0412, p-val = 0.8391

Supporting Tables

Table S3

TableS3 <- data.frame(question = rep(c("Q1", "Q2"), each =  2),
                      step = rep(c("Step 1", "Step 2"),2),
                      slope = c(round(d1s1$b[2],3), round(d1s2$b[2],3), round(d2s1$b[2],3), round(d2s2$b[2],3)),
                      se = c(round(d1s1$se[2],3), round(d1s2$se[2],3), round(d2s1$se[2],3), round(d2s2$se[2],3)),
                      pval =c(round(d1s1$pval[2],3), round(d1s2$pval[2],3), round(d2s1$pval[2],3), round(d2s2$pval[2],3)))

TableS3 %>%
  kbl() %>%
  kable_material(c("striped", "hover"), full_width = F)
question step slope se pval
Q1 Step 1 -1.546 0.206 0.000
Q1 Step 2 -0.871 0.099 0.000
Q2 Step 1 -0.544 0.165 0.001
Q2 Step 2 -0.366 0.089 0.000

Table S4

TableS4 <- data.frame(Stressor_type = Table2$Stressor_type,
                     Response_trait = Table2$Response_trait,
                     Overall_mean = Table2$Overall_mean,
                     Lower_95 = Table2$Lower_95,
                     Upper_95 = Table2$Upper_95,
                     P_value = Table2$P_value,
                     Overall_mean_adjusted = round(Q1.B$beta[-1],3),
                     Lower_95_adjusted = round(Q1.B$ci.lb[-1],3),
                     Upper_95_adjusted = round(Q1.B$ci.ub[-1],3),
                     P_value_adjusted = round(Q1.B$pval[-1],3),
                     N_effects = Table2$N_effects,
                     N_experiments = Table2$N_experiments,
                     N_host_taxa = Table2$N_host_taxa)

TableS4 %>%
  kbl() %>%
  kable_material(c("striped", "hover"), full_width = T)
Stressor_type Response_trait Overall_mean Lower_95 Upper_95 P_value Overall_mean_adjusted Lower_95_adjusted Upper_95_adjusted P_value_adjusted N_effects N_experiments N_host_taxa
Endogenous environment Fecundity -0.326 -0.802 0.150 0.180 -0.098 -0.555 0.359 0.675 14 3 3
Endogenous environment Survivorship -0.516 -0.715 -0.317 0.000 -0.176 -0.393 0.041 0.112 189 45 34
Chemical pollution Fecundity -0.141 -0.759 0.477 0.656 0.059 -0.553 0.671 0.850 14 4 2
Chemical pollution Survivorship -0.351 -0.640 -0.062 0.017 -0.106 -0.409 0.197 0.493 76 23 15
Resource limitation Fecundity -0.864 -1.287 -0.440 0.000 -0.694 -1.121 -0.267 0.001 59 6 4
Resource limitation Survivorship 0.002 -0.399 0.404 0.990 0.305 -0.093 0.703 0.134 32 6 5

Table S5

TableS5 <- data.frame(Stressor_type = Table3$Stressor_type,
                     Response_trait = Table3$Response_trait,
                     Overall_mean = Table3$Overall_mean,
                     Lower_95 = Table3$Lower_95,
                     Upper_95 = Table3$Upper_95,
                     P_value = Table3$P_value,
                     Overall_mean_adjusted = round(Q2.B$beta[-1],3),
                     Lower_95_adjusted = round(Q2.B$ci.lb[-1],3),
                     Upper_95_adjusted = round(Q2.B$ci.ub[-1],3),
                     P_value_adjusted = round(Q2.B$pval[-1],3),
                     N_effects = Table3$N_effects,
                     N_experiments = Table3$N_experiments,
                     N_host_taxa = Table3$N_host_taxa)

TableS5 %>%
  kbl() %>%
  kable_material(c("striped", "hover"), full_width = T)
Stressor_type Response_trait Overall_mean Lower_95 Upper_95 P_value Overall_mean_adjusted Lower_95_adjusted Upper_95_adjusted P_value_adjusted N_effects N_experiments N_host_taxa
Endogenous environment Prevalence 0.101 -0.152 0.354 0.435 0.255 -0.005 0.515 0.055 76 28 22
Endogenous environment Intensity 0.398 0.214 0.582 0.000 0.502 0.314 0.689 0.000 152 55 41
Endogenous environment Fecundity 0.032 -0.517 0.581 0.910 0.103 -0.437 0.642 0.710 9 4 4
Endogenous environment Survivorship -0.529 -0.705 -0.352 0.000 -0.405 -0.588 -0.222 0.000 162 66 47
Chemical pollution Prevalence -0.523 -1.017 -0.029 0.038 -0.395 -0.886 0.096 0.115 19 9 7
Chemical pollution Intensity -0.019 -0.267 0.229 0.879 0.050 -0.197 0.297 0.691 74 30 20
Chemical pollution Fecundity 0.089 -0.434 0.612 0.740 0.203 -0.315 0.722 0.442 16 7 4
Chemical pollution Survivorship -0.263 -0.521 -0.004 0.046 -0.151 -0.411 0.109 0.255 70 32 22
Resource limitation Prevalence -0.150 -0.633 0.333 0.542 -0.040 -0.518 0.438 0.870 15 5 4
Resource limitation Intensity -0.435 -0.792 -0.077 0.017 -0.331 -0.687 0.025 0.068 35 13 10
Resource limitation Fecundity -0.742 -1.138 -0.346 0.000 -0.676 -1.068 -0.285 0.001 31 8 5
Resource limitation Survivorship 0.005 -0.379 0.388 0.980 0.108 -0.273 0.489 0.578 27 10 9

R1

Here we conduct analyses suggested by reviewers.

Comment 1.R1 - Do infected and uninfected hosts differ in fitness?

We’ll address this question by Reviewer 1 in the control condition (i.e. in the absence of stressors)

Re arrange data

# columns to keep
keepers <- c( "ID", "Study", "Experiment", "Entered.by", "Host", "Host.type", "Parasite", "Parasite.type", "Study.type", "Trait", "Trait.type", "Statistics", "Level_C", "N_C", "Mean_C", "Variation_C", "Lwr_C", "Upr_C", "Success_C", "Variation.Type")

# data for uninfected hosts
dat_U <- dat[dat$Trait.category == "Uninfected demographic", keepers]

# change "Control" to "Uninfected"
colnames(dat_U) <- gsub("_C", "_U", colnames(dat_U))

# data for infected hosts
dat_I <- dat[dat$Trait.category == "Infected demographic", keepers]

# change "Control" to "Infected"
colnames(dat_I) <- gsub("_C", "_I", colnames(dat_I))

# merge data sets
dat_UI <- left_join(dat_U, dat_I, multiple = "all", by = c("Study", "Experiment", "Entered.by", "Host", "Host.type", "Parasite", "Parasite.type", "Study.type", "Trait.type", "Statistics", "Variation.Type"))

# remove experiments without infected data
dat_UI <- dat_UI[!is.na(dat_UI$Level_I),]

# how many effects
nrow(dat_UI)
## [1] 685

Calculate OR and approximate variance from contingency tables

N <- nrow(dat_UI)

for (i in 1:N) {
  # no correction need if none of the categories are 0
  if (dat_UI$Variation.Type[i] == "OR" &
      dat_UI$Success_U[i] > 0 &
      dat_UI$Success_I[i] > 0 &
      (dat_UI$N_U[i] - dat_UI$Success_U[i]) > 0 &
      (dat_UI$N_I[i] - dat_UI$Success_I[i]) > 0) {
    dat_UI$OR[i] <-
      (dat_UI$Success_I[i] * (dat_UI$N_U[i] - dat_UI$Success_U[i])) / (dat_UI$Success_U[i] * (dat_UI$N_I[i] - dat_UI$Success_I[i]))
    dat_UI$Log.OR[i] <- log(dat_UI$OR[i])
    # approximate variance
    dat_UI$Log.OR.v[i] <-
      1 / dat_UI$Success_I[i] + 1 / (dat_UI$N_I[i] - dat_UI$Success_I[i]) + 1 / dat_UI$Success_U[i] + 1 / (dat_UI$N_U[i] - dat_UI$Success_U[i])
  }
  else {
    # if at least one category is 0, we apply Yate's correction to avoid dividing by 0
    if (dat_UI$Variation.Type[i] == "OR" &
        (dat_UI$Success_U[i] == 0 |
        dat_UI$Success_I[i] == 0 |
        (dat_UI$N_U[i] - dat_UI$Success_U[i]) == 0 |
        (dat_UI$N_I[i] - dat_UI$Success_I[i]) == 0)) {
      dat_UI$OR[i] <-
        ((dat_UI$Success_I[i] + 0.5) * (dat_UI$N_U[i] - dat_UI$Success_U[i] + 0.5)) / ((dat_UI$Success_U[i] + 0.5) * (dat_UI$N_I[i] - dat_UI$Success_I[i] + 0.5))
      dat_UI$Log.OR[i] <- log(dat_UI$OR[i])
      # approximate variance
      dat_UI$Log.OR.v[i] <-
        1 / (dat_UI$Success_I[i] + 0.5)  + 1 / (dat_UI$N_I[i] - dat_UI$Success_I[i] + 0.5) + 1 / (dat_UI$Success_U[i] + 0.5) + 1 / (dat_UI$N_U[i] - dat_UI$Success_U[i] + 0.5)
    }
    else{
      # if there is no OR dat_UIa, make these columns NA
      dat_UI$OR[i] <- NA
      dat_UI$Log.OR[i] <- NA
      dat_UI$Log.OR.v[i] <- NA
    }
  }
}

Calculate SD for comparisons with continuous normally distributed variables

# create SD variables
dat_UI$SD_U <- vector(length = N)
dat_UI$SD_I <- vector(length = N)

for (i in 1:N) {
  # Calculate SD if SE is reported
  if (dat_UI$Variation.Type[i] == "SE") {
    dat_UI$SD_U[i] <- dat_UI$Variation_U[i] * sqrt(dat_UI$N_U[i])
    dat_UI$SD_I[i] <- dat_UI$Variation_I[i] * sqrt(dat_UI$N_I[i])
  } else {
    # calculate SD if a 95% confidence interval for a normal distribution is reported
    # this also applies to the Wald confidence interval for proportions, as it is a normal approximation to the binomial
    if (dat_UI$Variation.Type[i] == "CI" |
        dat_UI$Variation.Type[i] == "Wald CI") {
      # calculate SE from lower and upper limits
      temp_U1 <- (dat_UI$Upr_U[i] - dat_UI$Mean_U[i]) / 1.96
      temp_U2 <- (dat_UI$Mean_U[i] - dat_UI$Lwr_U[i]) / 1.96
      
      temp_I1 <- (dat_UI$Upr_I[i] - dat_UI$Mean_I[i]) / 1.96
      temp_I2 <- (dat_UI$Mean_I[i] - dat_UI$Lwr_I[i]) / 1.96
      
      # average digitized/recorded values of SE (because we have two) and transform to SD
      dat_UI$SD_U[i] <-
        mean(abs(c(temp_U1, temp_U2))) * sqrt(dat_UI$N_U[i])
      dat_UI$SD_I[i] <-
        mean(abs(c(temp_U1, temp_U2))) * sqrt(dat_UI$N_I[i])
    } else {
      # approximate SD if IQ range is reported
      if (dat_UI$Variation.Type[i] == "IQ") {
        dat_UI$SD_U[i] <- (dat_UI$Upr_U[i] - dat_UI$Lwr_U[i]) / 1.35
        dat_UI$SD_I[i] <- (dat_UI$Upr_I[i] - dat_UI$Lwr_I[i]) / 1.35
      }
      else {
        # if SD is reported leave as such
        if (dat_UI$Variation.Type[i] == "SD") {
          dat_UI$SD_U[i] <- dat_UI$Variation_U[i]
          dat_UI$SD_I[i] <- dat_UI$Variation_I[i]
        } else {
          # if there is no appropriate dat_UIa, make these columns NA
          dat_UI$SD_U[i] <- NA
          dat_UI$SD_I[i] <- NA
        }
      }
    }
  }
}

Get effect sizes and variances

# within groups standard deviation
dat_UI$S_within <-
  sqrt(((dat_UI$N_U - 1) * dat_UI$SD_U ^ 2 + (dat_UI$N_I - 1) * dat_UI$SD_I ^ 2) / (dat_UI$N_U + dat_UI$N_I - 2))

# standardized effect size
for (i in 1:N) {
  # standardized mean difference
  # we can't include data without variance
  if (is.na(dat_UI$S_within[i]) == FALSE & dat_UI$S_within[i] > 0) {
    dat_UI$d[i] <- (dat_UI$Mean_I[i] - dat_UI$Mean_U[i]) / dat_UI$S_within[i]
    dat_UI$var.d[i] <-
      (dat_UI$N_U[i] + dat_UI$N_I[i]) / (dat_UI$N_U[i] * dat_UI$N_I[i]) + (dat_UI$d[i] ^ 2) /
      (2 * (dat_UI$N_U[i] + dat_UI$N_I[i]))
  }
  else{
    # if only F statistic is reported - no cases of this yet!
    if (dat_UI$Variation.Type[i] == "F_I") {
      dat_UI$d[i] <-
        sqrt(dat_UI$F_I[i] * (dat_UI$N_I[i] + dat_UI$N_U[i]) / (dat_UI$N_U[i] * dat_UI$N_I[i]))
      dat_UI$var.d[i] <-
        (dat_UI$N_U[i] + dat_UI$N_I[i]) / (dat_UI$N_U[i] * dat_UI$N_I[i]) + (dat_UI$d[i] ^ 2) /
        (2 * (dat_UI$N_U[i] + dat_UI$N_I[i]))
      
    } else{
      # if Z and N are reported for a regression analysis
      if (dat_UI$Variation.Type[i] == "Z_reg") {
        r <- dat_UI$Z_I[i] / sqrt(dat_UI$Z_N[i])
        dat_UI$d[i] <- (2 * r) / (sqrt(1 - r ^ 2))
        Vr <- (1 - r ^ 2) ^ 2 / (dat_UI$Z_N[i] - 1)
        dat_UI$var.d[i] <- (4 * Vr) / (1 - r ^ 2) ^ 3
      }
      else{
        # if t and N are reported for a regression analysis
        if (dat_UI$Variation.Type[i] == "t") {
          r <- sqrt((dat_UI$Z_I[i] ^ 2) / (dat_UI$Z_I[i] ^ 2 + dat_UI$Z_df[i]))
          dat_UI$d[i] <- (2 * r) / (sqrt(1 - r ^ 2))
          Vr <- (1 - r ^ 2) ^ 2 / (dat_UI$R_N[i] - 1)
          dat_UI$var.d[i] <- (4 * Vr) / (1 - r ^ 2) ^ 3
        }
        else{
          # If Z is reported from a comparison of two independent groups
          if (dat_UI$Variation.Type[i] == "Z") {
            dat_UI$d[i] <-
              sqrt(abs(dat_UI$Z_I[i]) * sqrt(dat_UI$N_U[i] + dat_UI$N_I[i]) / (1 - sqrt(
                dat_UI$Z_I[i] ^ 2 * (dat_UI$N_U[i] + dat_UI$N_I[i]) ^ -1
              )))
            dat_UI$var.d[i] <-
              (dat_UI$N_U[i] + dat_UI$N_I[i]) / (dat_UI$N_U[i] * dat_UI$N_I[i]) + (dat_UI$d[i] ^ 2) /
              (2 * (dat_UI$N_U[i] + dat_UI$N_I[i]))
          }
          else{
            # converting log OR
            if (is.na(dat_UI$Log.OR.v[i]) == FALSE &
                dat_UI$Log.OR.v[i] > 0) {
              dat_UI$d[i] <- dat_UI$Log.OR[i] * (sqrt(3) / pi)
              dat_UI$var.d[i] <- dat_UI$Log.OR.v[i] * 3 / pi ^ 2
            } else {
              # for now, leave other types of effects as NA
              dat_UI$d[i] <- NA
              dat_UI$var.d[i] <- NA
            }
          }
        }
      }
    }
  }
}

Correct for sample size and get g

# sample size correction factor
for (i in 1:N) {
  # if odds ratio or comparison between means
  if (is.na(dat_UI$N_U[i]) == F) {
    dat_UI$J[i] <- 1 - 3 / (4 * (dat_UI$N_U[i] + dat_UI$N_I[i] - 2) - 1)
  } else {
    # if correlation
    if (is.na(dat_UI$R_N[i]) == F) {
      dat_UI$J[i] <- 1 - 3 / (4 * (dat_UI$R_N[i] - 2) - 1)
    }else {
      dat_UI$J[i] <- NA
    }
  }
}
# corrected effect size
dat_UI$g <- dat_UI$J * dat_UI$d

# and variance
dat_UI$var.g <- (dat_UI$J ^ 2) * dat_UI$var.d

# remove studies without variance for which g had been set to NA
dat_UI <- dat_UI[!is.na(dat_UI$g),]

# remove duplicated rows which were created whenever controls were used more than once in main analysis.
dat_UI <- dat_UI[!duplicated(dat_UI[,-c(1,21)]), ]

Flip effect sizes for remaining mortality effects

mort <- grep(pattern = "[M,m]ort", x = dat_UI$Trait.x)

for(i in mort){
  dat_UI$g[i] <- -dat_UI$g[i]
}

Identify negative eigenvalues

As in our main analysis, we find out which studies have negative eigenvalues in their variance-covar matrix. This means covariance between treatment effects with the same control is larger than variance and it’s an indicator of suspiciously low variance and potentially an incorrectly reported N or an error while digitizing data.

studies_1 <- unique(dat_UI$Study)
studies_to_Uheck <- c()
negative_eigen <- c()

for (k in studies_1) {
  dat_study <- dat_UI[which(dat_UI$Study == k),]
  varcovmat = matrix(0,
                     nrow = dim(dat_study)[1],
                     ncol = dim(dat_study)[1])
  
  for (i in 1:dim(dat_study)[1]) {
    for (j in 1:dim(dat_study)[1]) {
      if (i == j) {
        varcovmat[i, j] = dat_study$var.g[i]
      } else{
        if (dat_study[i, "Experiment"] == dat_study[j, "Experiment"] &
            dat_study[i, "Level_U"] == dat_study[j, "Level_U"]) {
          varcovmat[i, j] = 1 / dat_study[i, "N_U"] + dat_study$g[i] * dat_study$g[j] /
            (dat_study[i, "N_U"] + dat_study[i, "N_I"] + dat_study[j, "N_I"])
        }
      }
    }
  }
  
  val <- eigen(varcovmat)
  for (m in 1:length(val$values)) {
    if (val$values[m] < 0) {
    studies_to_Uheck <- append(studies_to_Uheck, k)
    negative_eigen <- append(negative_eigen, val$values[m])
    }
  }
}
    
Neg_eigen_UI <- data.frame(study = studies_to_Uheck, eigen = negative_eigen)

Neg_eigen_UI %>%
  kbl() %>%
  kable_material(c("striped", "hover"), full_width = F)
study eigen
Ashraf et al. 2017 Environ Sci Pollut Res -0.0181032
Ashraf et al. 2017 Environ Sci Pollut Res -0.1759683
Ashraf et al. 2017 Environ Sci Pollut Res -0.3098886
Ashraf et al. 2017 Environ Sci Pollut Res -0.8423039
Manzi et al. 2019. Freshwater Biology -0.0103166
Manzi et al. 2019. Freshwater Biology -0.0327272
# exclude suspicious studies with negative eigenvalues and which already raised flags in the main analysis
dat_UI <- dat_UI[-grep("Ashraf et al. 2017 Environ Sci Pollut Res", dat_UI$Study),]
dat_UI <- dat_UI[-grep("Civitello et al. 2020. Proc B", dat_UI$Study),]

# I'm also excluding Hall et al. as they also had too negative eigen values
dat_UI <- dat_UI[-grep("Hall et al. 2013. Oecologia", dat_UI$Study),]
dat_UI <- dat_UI[-grep("Manzi et al. 2019. Freshwater Biology", dat_UI$Study),]
        
    
varcovmat_UI = matrix(0, nrow = dim(dat_UI)[1], ncol = dim(dat_UI)[1])

for (i in 1:dim(dat_UI)[1]) {
  for (j in 1:dim(dat_UI)[1]) {
    if (i == j) {varcovmat_UI[i,j] = dat_UI$var.g[i]}else{
      if (dat_UI[i, "Experiment"] == dat_UI[j, "Experiment"] & dat_UI[i, "Level_U"] == dat_UI[j, "Level_U"]) {
        varcovmat_UI[i,j] = 1/dat_UI[i,"N_U"] + dat_UI$g[i]*dat_UI$g[j]/(dat_UI[i,"N_U"] + dat_UI[i,"N_I"] + dat_UI[j,"N_I"]) 
      }
    }
  }
}

# correct negative eigens in the few studies with large covar to var ratios
varcovmat_UI_PD <- nearPD(varcovmat_UI,  keepDiag = TRUE)

Model without moderators

We fit the model to test whether infection reduces fitness (fecundity or survivorship)

Q0.R1.f <-
  rma.mv(
    g ~ Trait.type -1,
    V = varcovmat_UI_PD$mat,
    random = list(~ 1 | Experiment,  ~ 1 | ID.x),
    data = dat_UI,
    method = "REML"
  )

pal <- RColorBrewer::brewer.pal(6,"Spectral")

Q0.R1.f_OP <-
  orchard_plot(
    Q0.R1.f,
    xlab = "Standardised mean difference (Hedge's g)",
    transfm = "none",
    angle = 0,
    alpha = 0.7
  ) +
  scale_fill_manual(values = rep(pal, each = 2)) + scale_colour_manual(values = rep(pal, each = 2)) +
  theme(legend.position = c(0.25, 0.01), axis.title = element_text(size = 15), axis.text.y = element_text(size = 15)) +
  scale_y_discrete(
    labels = c(
      "Fecundity",
      "Survivorship"
    )
  )

Q0.R1.f_OP

Forest plot

summary(Q0.R1.f)
## 
## Multivariate Meta-Analysis Model (k = 308; method: REML)
## 
##     logLik    Deviance         AIC         BIC        AICc  ​ 
## -1161.8670   2323.7340   2331.7340   2346.6284   2331.8669   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed      factor 
## sigma^2.1  0.9933  0.9966     72     no  Experiment 
## sigma^2.2  0.0259  0.1608    118     no        ID.x 
## 
## Test for Residual Heterogeneity:
## QE(df = 306) = 3629.6611, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## QM(df = 2) = 56.5456, p-val < .0001
## 
## Model Results:
## 
##                         estimate      se     zval    pval    ci.lb    ci.ub​ 
## Trait.typeFecundity      -2.0056  0.2984  -6.7215  <.0001  -2.5905  -1.4208 
## Trait.typeSurvivorship   -0.6775  0.1428  -4.7457  <.0001  -0.9573  -0.3977 
##  
## Trait.typeFecundity     *** 
## Trait.typeSurvivorship  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest.default(x= Q0.R1.f$beta, sei =  Q0.R1.f$se, ci.lb =  Q0.R1.f$ci.lb, ci.ub =  Q0.R1.f$ci.ub, 
                annotate=TRUE, showweights=T, header=F,
                slab = c("Fecundity", "Survivorship"))

Comment 2.R1 - Do aquatic and terrestrial systems differ in their responses to stressors?

To address this question by Reviewer 2, we’ll add an ‘Environment’ moderator to the final models in Q1 and Q2, presented in the main text. The results of these models are presented in the Supporting Material.

Q1: Best Model + Environment

Effects in survival and fecundity under different types of stress and in aquatic vs terrestrial habitats, combining infected and uninfected.

Q1.R1 <-
  rma.mv(
    g ~ Environment:Trait.type:Gradient.category -1 ,
    V = varcovmat_1_PD$mat,
    random = list( ~ 1 | Experiment, ~ 1 | ID),
    data = dat_1,
    method = "REML"
  )
## Warning: Redundant predictors dropped from the model.
summary(Q1.R1)
## 
## Multivariate Meta-Analysis Model (k = 384; method: REML)
## 
##    logLik   Deviance        AIC        BIC       AICc  ​ 
## -540.8105  1081.6210  1103.6210  1146.8172  1104.3483   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed      factor 
## sigma^2.1  0.2203  0.4693     78     no  Experiment 
## sigma^2.2  0.3630  0.6025    384     no          ID 
## 
## Test for Residual Heterogeneity:
## QE(df = 375) = 214780.6070, p-val < .0001
## 
## Test of Moderators (coefficients 1:9):
## QM(df = 9) = 57.7444, p-val < .0001
## 
## Model Results:
## 
##                                                                             estimate​ 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment          -0.3890 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment       -0.6141 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment    0.1304 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution            -0.0724 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution         -0.1986 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution     -0.4953 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource             -0.9815 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource         -0.4179 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource          -0.0528 
##                                                                                 se 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment         0.2419 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment      0.1048 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment  0.2793 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution           0.3140 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution        0.1982 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution    0.2038 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource            0.2331 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource        0.4802 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource         0.2039 
##                                                                                zval 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment         -1.6084 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment      -5.8582 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment   0.4668 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution           -0.2306 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution        -1.0019 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution    -2.4305 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource            -4.2102 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource        -0.8702 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource         -0.2588 
##                                                                               pval 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment         0.1078 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment      <.0001 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment  0.6407 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution           0.8176 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution        0.3164 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution    0.0151 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource            <.0001 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource        0.3842 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource         0.7958 
##                                                                               ci.lb 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment         -0.8630 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment      -0.8196 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment  -0.4170 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution           -0.6879 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution        -0.5871 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution    -0.8947 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource            -1.4384 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource        -1.3591 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource         -0.4523 
##                                                                               ci.ub 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment          0.0850 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment      -0.4086 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment   0.6777 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution            0.5430 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution         0.1899 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution    -0.0959 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource            -0.5246 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource         0.5233 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource          0.3468 
##  
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment      *** 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution      * 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource            *** 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest.default(
  x = Q1.R1$beta,
  sei =  Q1.R1$se,
  ci.lb =  Q1.R1$ci.lb,
  ci.ub =  Q1.R1$ci.ub,
  annotate = TRUE,
  showweights = T,
  header = F,
  slab = c(
    "Aquatic environment fecundity",
    "Aquatic environment survivorship",
    "Terrestrial environment survivorship",
    "Aquatic pollution  fecundity",
    "Aquatic pollution survivorship",
    "Terrestrial pollution  survivorship",
    "Aquatic resource  fecundity",
    "Terrestrial resource  fecundity",
    "Aquatic resource survivorship"
  )
)

pal <- RColorBrewer::brewer.pal(11,"Spectral")

Q1.R1_OP <-
  orchard_plot(
    Q1.R1,
    xlab = "Standardised mean difference (Hedge's g)",
    transfm = "none",
    angle = 0,
    alpha = 0.7
  ) +
  scale_fill_manual(values = pal[c(7,8,8,10,11,11,5,5,4)]) +
  scale_colour_manual(values = pal[c(7,8,8,10,11,11,5,5,4)]) +
  theme(legend.position = c(0.25, 0.01), axis.title = element_text(size = 15), axis.text.y = element_text(size = 15)) +
  scale_y_discrete(
    labels = c("AQ:EE:fecundity",
               "AQ:EE:survivorship",
               "TE:EE:survivorship",
               "AQ:CP:fecundity",
               "AQ:CP:survivorship",
               "TE:CP:survivorship",
               "AQ:RL:fecundity",
               "TE:RL:fecundity",
               "AQ:RL:survivorship"))
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
Q1.R1_OP

Make Table with results from Q1.R1

TableQ1.R1 <-
  data.frame(
    Environment = c(
      "Aquatic",
      "Aquatic",
      "Terrestrial",
      "Aquatic",
      "Aquatic",
      "Terrestrial",
      "Aquatic",
      "Terrestrial",
      "Aquatic"
    ),
    Gradient.category = rep(
      c(
        "Environment",
        "Pollution",
        "Resource"
      ),
      each = 3
    ),
    Trait.type = c(
      "Fecundity",
      "Survivorship",
      "Survivorship",
      "Fecundity",
      "Survivorship",
      "Survivorship",
      "Fecundity",
      "Fecundity",
      "Survivorship"
    ),
    Overall_mean = round(Q1.R1$beta, 3),
    Lower_95 = round(Q1.R1$ci.lb, 3),
    Upper_95 = round(Q1.R1$ci.ub, 3),
    P_value = round(Q1.R1$pval, 3))

N_effects <- dat_1 %>% group_by(Environment, Gradient.category, Trait.type) %>% summarise(N_effects = length(ID))
## `summarise()` has grouped output by 'Environment', 'Gradient.category'. You can
## override using the `.groups` argument.
N_experiments <- dat_1 %>% group_by(Environment, Gradient.category, Trait.type) %>% summarise(N_experiments = length(unique(Experiment)))
## `summarise()` has grouped output by 'Environment', 'Gradient.category'. You can
## override using the `.groups` argument.
N_hosts <- dat_1 %>% group_by(Environment, Gradient.category, Trait.type) %>% summarise(N_hosts = length(unique(Host)))
## `summarise()` has grouped output by 'Environment', 'Gradient.category'. You can
## override using the `.groups` argument.
TableQ1.R1 <-
  left_join(TableQ1.R1, N_effects) %>% left_join(N_experiments) %>%   left_join(N_experiments) %>%
  dplyr::rename(Stressor_type = Gradient.category,
                Response_trait = Trait.type) %>%
  mutate(Stressor_type = plyr::revalue(
    Stressor_type,
    c(
      "Environment" = "Endogenous Environment",
      "Pollution" = "Chemical Pollution",
      "Resource" = "Resource Limitation"
    )
  ))
## Joining with `by = join_by(Environment, Gradient.category, Trait.type)`
## Joining with `by = join_by(Environment, Gradient.category, Trait.type)`
## Joining with `by = join_by(Environment, Gradient.category, Trait.type,
## N_experiments)`
rownames(TableQ1.R1) <- NULL

TableQ1.R1 %>%
  kbl() %>%
  kable_material(c("striped", "hover"), full_width = F)
Environment Stressor_type Response_trait Overall_mean Lower_95 Upper_95 P_value N_effects N_experiments
Aquatic Endogenous Environment Fecundity -0.389 -0.863 0.085 0.108 14 3
Aquatic Endogenous Environment Survivorship -0.614 -0.820 -0.409 0.000 170 40
Terrestrial Endogenous Environment Survivorship 0.130 -0.417 0.678 0.641 19 5
Aquatic Chemical Pollution Fecundity -0.072 -0.688 0.543 0.818 14 4
Aquatic Chemical Pollution Survivorship -0.199 -0.587 0.190 0.316 44 12
Terrestrial Chemical Pollution Survivorship -0.495 -0.895 -0.096 0.015 32 11
Aquatic Resource Limitation Fecundity -0.981 -1.438 -0.525 0.000 23 5
Terrestrial Resource Limitation Fecundity -0.418 -1.359 0.523 0.384 36 1
Aquatic Resource Limitation Survivorship -0.053 -0.452 0.347 0.796 32 6

Does accounting for the transmission environment reduce heterogeneity? We calculate I2 using the formula above from Nakagawa and Santos (2012) for the model without the var-covar matrix.

Q1.R1.R <-
  rma.mv(
    g ~ Environment:Trait.type:Gradient.category -1 ,
    V = dat_1$var.g,
    random = list( ~ 1 | Experiment, ~ 1 | ID),
    data = dat_1,
    method = "REML"
  )
## Warning: Redundant predictors dropped from the model.
I2.Q1.R1 <- I2(Q1.R1.R, v = dat_1$var.g, obs = "ID")*100
## Warning in class(model) != "rma.mv" && class(model) != "rma": 'length(x) = 2 >
## 1' in coercion to 'logical(1)'
I2.Q1.R1
##            I2_Est. 2.5% CI 97.5% CI
## Experiment   51.52   42.40    59.52
## total        93.63   92.42    94.63

Q2: Best Model + Environment

Trait type x Gradient category x Transmission environment interaction

Q2.R1 <-
  rma.mv(
    g ~ Environment:Trait.type:Gradient.category -1,
    V = varcovmat_2_PD$mat,
    random = list(~ 1 | Experiment,  ~ 1 | ID),
    data = dat_2,
    method = "REML",
)

summary(Q2.R1)
## 
## Multivariate Meta-Analysis Model (k = 686; method: REML)
## 
##    logLik   Deviance        AIC        BIC       AICc  ​ 
## -999.5082  1999.0163  2051.0163  2167.8932  2053.2273   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed      factor 
## sigma^2.1  0.1790  0.4231    113     no  Experiment 
## sigma^2.2  0.4244  0.6515    686     no          ID 
## 
## Test for Residual Heterogeneity:
## QE(df = 662) = 1419271.1848, p-val < .0001
## 
## Test of Moderators (coefficients 1:24):
## QM(df = 24) = 147.1773, p-val < .0001
## 
## Model Results:
## 
##                                                                             estimate​ 
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryEnvironment          0.3076 
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryEnvironment      0.0876 
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryEnvironment           0.4751 
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryEnvironment      -0.0980 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment           0.0299 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryEnvironment       0.1172 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment       -0.7433 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment    0.1154 
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryPollution           -1.0873 
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryPollution       -0.2194 
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryPollution             0.0173 
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryPollution        -0.0897 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution            -0.1289 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryPollution         0.2036 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution         -0.2680 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution     -0.2483 
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryResource            -0.1147 
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryResource        -0.8293 
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryResource             -0.5378 
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryResource         -0.1989 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource             -0.7288 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource         -0.7240 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource           0.0670 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryResource      -0.2270 
##                                                                                 se 
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryEnvironment        0.1629 
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryEnvironment    0.2250 
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryEnvironment         0.1009 
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryEnvironment     0.2509 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment         0.3119 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryEnvironment     0.5969 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment      0.1019 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment  0.1986 
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryPollution          0.4160 
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryPollution      0.3152 
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryPollution           0.1640 
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryPollution       0.1999 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution           0.3957 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryPollution       0.3577 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution        0.1771 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution    0.1963 
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryResource           0.2557 
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryResource       0.9515 
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryResource            0.2155 
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryResource        0.3399 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource            0.2431 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource        0.3619 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource         0.2220 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryResource     0.4022 
##                                                                                zval 
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryEnvironment         1.8881 
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryEnvironment     0.3895 
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryEnvironment          4.7067 
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryEnvironment     -0.3907 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment          0.0959 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryEnvironment      0.1964 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment      -7.2947 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment   0.5809 
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryPollution          -2.6135 
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryPollution      -0.6962 
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryPollution            0.1053 
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryPollution       -0.4488 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution           -0.3259 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryPollution        0.5691 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution        -1.5130 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution    -1.2647 
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryResource           -0.4487 
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryResource       -0.8716 
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryResource            -2.4961 
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryResource        -0.5853 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource            -2.9984 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource        -2.0003 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource          0.3019 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryResource     -0.5644 
##                                                                               pval 
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryEnvironment        0.0590 
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryEnvironment    0.6969 
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryEnvironment         <.0001 
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryEnvironment     0.6960 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment         0.9236 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryEnvironment     0.8443 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment      <.0001 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment  0.5613 
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryPollution          0.0090 
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryPollution      0.4863 
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryPollution           0.9161 
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryPollution       0.6536 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution           0.7445 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryPollution       0.5693 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution        0.1303 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution    0.2060 
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryResource           0.6537 
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryResource       0.3834 
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryResource            0.0126 
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryResource        0.5583 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource            0.0027 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource        0.0455 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource         0.7627 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryResource     0.5725 
##                                                                               ci.lb 
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryEnvironment        -0.0117 
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryEnvironment    -0.3534 
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryEnvironment          0.2773 
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryEnvironment     -0.5898 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment         -0.5815 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryEnvironment     -1.0526 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment      -0.9430 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment  -0.2739 
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryPollution          -1.9027 
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryPollution      -0.8371 
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryPollution           -0.3041 
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryPollution       -0.4814 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution           -0.9045 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryPollution       -0.4975 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution        -0.6151 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution    -0.6331 
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryResource           -0.6159 
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryResource       -2.6941 
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryResource            -0.9601 
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryResource        -0.8651 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource            -1.2052 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource        -1.4334 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource         -0.3680 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryResource     -1.0152 
##                                                                               ci.ub 
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryEnvironment         0.6268 
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryEnvironment     0.5286 
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryEnvironment          0.6730 
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryEnvironment      0.3937 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment          0.6412 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryEnvironment      1.2871 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment      -0.5436 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment   0.5047 
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryPollution          -0.2719 
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryPollution       0.3983 
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryPollution            0.3386 
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryPollution        0.3020 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution            0.6466 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryPollution        0.9046 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution         0.0792 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution     0.1365 
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryResource            0.3864 
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryResource        1.0356 
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryResource            -0.1155 
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryResource         0.4672 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource            -0.2524 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource        -0.0146 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource          0.5020 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryResource      0.5613 
##  
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryEnvironment          . 
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryEnvironment 
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryEnvironment         *** 
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryEnvironment 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryEnvironment 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryEnvironment 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryEnvironment      *** 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryEnvironment 
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryPollution           ** 
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryPollution 
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryPollution 
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryPollution 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryPollution 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryPollution 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryPollution 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryPollution 
## EnvironmentAquatic:Trait.typePrevalence:Gradient.categoryResource 
## EnvironmentTerrestrial:Trait.typePrevalence:Gradient.categoryResource 
## EnvironmentAquatic:Trait.typeIntensity:Gradient.categoryResource              * 
## EnvironmentTerrestrial:Trait.typeIntensity:Gradient.categoryResource 
## EnvironmentAquatic:Trait.typeFecundity:Gradient.categoryResource             ** 
## EnvironmentTerrestrial:Trait.typeFecundity:Gradient.categoryResource          * 
## EnvironmentAquatic:Trait.typeSurvivorship:Gradient.categoryResource 
## EnvironmentTerrestrial:Trait.typeSurvivorship:Gradient.categoryResource 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest.default(
  x = Q2.R1$beta,
  sei =  Q2.R1$se,
  ci.lb =  Q2.R1$ci.lb,
  ci.ub =  Q2.R1$ci.ub,
  annotate = TRUE,
  showweights = T,
  header = F,
  slab = c(
    "AQ:EE:prevalence",
    "TE:EE:prevalence",
    "AQ:EE:intensity",
    "TE:EE:intensity",
    "AQ:EE:fecundity",
    "TE:EE:fecundity",
    "AQ:EE:survivorship",
    "TE:EE:survivorship",
    "AQ:CP:prevalence",
    "TE:CP:prevalence",
    "AQ:CP:intensity",
    "TE:CP:intensity",
    "AQ:CP:fecundity",
    "TE:CP:fecundity",
    "AQ:CP:survivorship",
    "TE:CP:survivorship",
    "AQ:RL:prevalence",
    "TE:RL:prevalence",
    "AQ:RL:intensity",
    "TE:RL:intensity",
    "AQ:RL:fecundity",
    "TE:RL:fecundity",
    "AQ:RL:survivorship",
    "TE:RL:survivorship"
  )
)

par_short <- c("AQ:EE:prevalence",
    "TE:EE:prevalence",
    "AQ:EE:intensity",
    "TE:EE:intensity",
    "AQ:EE:fecundity",
    "TE:EE:fecundity",
    "AQ:EE:survivorship",
    "TE:EE:survivorship",
    "AQ:CP:prevalence",
    "TE:CP:prevalence",
    "AQ:CP:intensity",
    "TE:CP:intensity",
    "AQ:CP:fecundity",
    "TE:CP:fecundity",
    "AQ:CP:survivorship",
    "TE:CP:survivorship",
    "AQ:RL:prevalence",
    "TE:RL:prevalence",
    "AQ:RL:intensity",
    "TE:RL:intensity",
    "AQ:RL:fecundity",
    "TE:RL:fecundity",
    "AQ:RL:survivorship",
    "TE:RL:survivorship")

pal <- RColorBrewer::brewer.pal(11,"Spectral")

Q2.R1_OP <-
  orchard_plot(
    Q2.R1,
    xlab = "Standardised mean difference (Hedge's g)",
    transfm = "none",
    angle = 0,
    alpha = 0.7
  ) +
  theme(legend.position = c(0.3, 0.01), axis.title = element_text(size = 15), axis.text.y = element_text(size = 15)) +
  scale_fill_manual(values = rep(pal[c(7,8,10,11,5,4)], each = 4)) +
  scale_colour_manual(values = rep(pal[c(7,8,10,11,5,4)], each = 4)) +
  scale_y_discrete(labels = par_short)
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
Q2.R1_OP

Make Table with results from Q2.R1

TableQ2.R1 <-
  data.frame(
    Environment = rep(c(
      "Aquatic",
      "Terrestrial"), 12),
    Gradient.category = rep(
      c(
        "Environment",
        "Pollution",
        "Resource"
      ),
      each = 8
    ),
    Trait.type = rep(c(
      "Prevalence",
      "Prevalence",
      "Intensity",
      "Intensity",
      "Fecundity",
      "Fecundity",
      "Survivorship",
      "Survivorship"), 3
    ),
    Overall_mean = round(Q2.R1$beta, 3),
    Lower_95 = round(Q2.R1$ci.lb, 3),
    Upper_95 = round(Q2.R1$ci.ub, 3),
    P_value = round(Q2.R1$pval, 3))

N_effects <- dat_2 %>% group_by(Environment, Gradient.category, Trait.type) %>% summarise(N_effects = length(ID))
## `summarise()` has grouped output by 'Environment', 'Gradient.category'. You can
## override using the `.groups` argument.
N_experiments <- dat_2 %>% group_by(Environment, Gradient.category, Trait.type) %>% summarise(N_experiments = length(unique(Experiment)))
## `summarise()` has grouped output by 'Environment', 'Gradient.category'. You can
## override using the `.groups` argument.
N_hosts <- dat_2 %>% group_by(Environment, Gradient.category, Trait.type) %>% summarise(N_hosts = length(unique(Host)))
## `summarise()` has grouped output by 'Environment', 'Gradient.category'. You can
## override using the `.groups` argument.
TableQ2.R1 <-
  left_join(TableQ2.R1, N_effects) %>% left_join(N_experiments) %>%   left_join(N_experiments) %>%
  dplyr::rename(Stressor_type = Gradient.category,
                Response_trait = Trait.type) %>%
  mutate(Stressor_type = plyr::revalue(
    Stressor_type,
    c(
      "Environment" = "Endogenous Environment",
      "Pollution" = "Chemical Pollution",
      "Resource" = "Resource Limitation"
    )
  ))
## Joining with `by = join_by(Environment, Gradient.category, Trait.type)`
## Joining with `by = join_by(Environment, Gradient.category, Trait.type)`
## Joining with `by = join_by(Environment, Gradient.category, Trait.type,
## N_experiments)`
rownames(TableQ2.R1) <- NULL

TableQ2.R1 %>%
  kbl() %>%
  kable_material(c("striped", "hover"), full_width = F)
Environment Stressor_type Response_trait Overall_mean Lower_95 Upper_95 P_value N_effects N_experiments
Aquatic Endogenous Environment Prevalence 0.308 -0.012 0.627 0.059 41 21
Terrestrial Endogenous Environment Prevalence 0.088 -0.353 0.529 0.697 35 7
Aquatic Endogenous Environment Intensity 0.475 0.277 0.673 0.000 135 50
Terrestrial Endogenous Environment Intensity -0.098 -0.590 0.394 0.696 17 5
Aquatic Endogenous Environment Fecundity 0.030 -0.581 0.641 0.924 7 3
Terrestrial Endogenous Environment Fecundity 0.117 -1.053 1.287 0.844 2 1
Aquatic Endogenous Environment Survivorship -0.743 -0.943 -0.544 0.000 125 56
Terrestrial Endogenous Environment Survivorship 0.115 -0.274 0.505 0.561 37 10
Aquatic Chemical Pollution Prevalence -1.087 -1.903 -0.272 0.009 8 4
Terrestrial Chemical Pollution Prevalence -0.219 -0.837 0.398 0.486 11 5
Aquatic Chemical Pollution Intensity 0.017 -0.304 0.339 0.916 48 17
Terrestrial Chemical Pollution Intensity -0.090 -0.481 0.302 0.654 26 13
Aquatic Chemical Pollution Fecundity -0.129 -0.905 0.647 0.745 7 4
Terrestrial Chemical Pollution Fecundity 0.204 -0.498 0.905 0.569 9 3
Aquatic Chemical Pollution Survivorship -0.268 -0.615 0.079 0.130 46 18
Terrestrial Chemical Pollution Survivorship -0.248 -0.633 0.136 0.206 24 14
Aquatic Resource Limitation Prevalence -0.115 -0.616 0.386 0.654 14 4
Terrestrial Resource Limitation Prevalence -0.829 -2.694 1.036 0.383 1 1
Aquatic Resource Limitation Intensity -0.538 -0.960 -0.116 0.013 25 9
Terrestrial Resource Limitation Intensity -0.199 -0.865 0.467 0.558 10 4
Aquatic Resource Limitation Fecundity -0.729 -1.205 -0.252 0.003 18 6
Terrestrial Resource Limitation Fecundity -0.724 -1.433 -0.015 0.045 13 2
Aquatic Resource Limitation Survivorship 0.067 -0.368 0.502 0.763 21 7
Terrestrial Resource Limitation Survivorship -0.227 -1.015 0.561 0.573 6 3

Does accounting for the transmission environment reduce heterogeneity? We calculate I2 using the formula above from Nakagawa and Santos (2012) for the model without the var-covar matrix.

Q2.R1.R <-
  rma.mv(
    g ~ Environment:Trait.type:Gradient.category -1 ,
    V = dat_2$var.g,
    random = list( ~ 1 | Experiment, ~ 1 | ID),
    data = dat_2,
    method = "REML"
  )
I2.Q2.R1 <- I2(Q2.R1.R, v = dat_1$var.g, obs = "ID")*100
## Warning in class(model) != "rma.mv" && class(model) != "rma": 'length(x) = 2 >
## 1' in coercion to 'logical(1)'
I2.Q2.R1
##            I2_Est. 2.5% CI 97.5% CI
## Experiment   29.89   24.21    35.83
## total        94.35   93.74    94.89